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Principal  Investigator:  Jin  Au  Kong 


FINAL  REPORT 


Under  the  sponsorship  of  the  ONR  Contract  N00014-86-K-0533  we  have  published 
15  referenced  journal  and  conference  papers  on  time  domain  electromagnetic  waves  in 
multilayer  media. '  ) 

^  A  rigorous  dyadic  Green's  function  formulation  in  the  spectral  domain  is  used  to 
study  the  dispersion  characteristics  of  signal  strip  lines  in  the  presence  of  metallic  crossing 
strips. 

V  A  set  of  coupled  vector  integral  equations  for  the  current  distribution  on  the  conduc¬ 
tors  is  derived.  Galerkin’s  method  is  then  applied  to  derive  the  matrix  eigenvalue  equation 
for  the  propagation  constant.  The  dispersion  properties  of  the  signal  lines  are  studied  for 
both  cases  of  finite  and  infinite  length  crossing  strips. 

^  The  effects  of  the  structure  dimensions  on  the  passband  and  stopband  characteristics 
are  investigated.  For  crossing  strips  of  finite  length,  the  stopband  is  mainly  affected  by  the 
period,  the  crossing  strip  length,  and  the  separation  between  the  signal  and  the  crossing 
strips.  For  crossing  strips  of  infinite  length  carrying  travelling  waves,  attenuation  along 
the  signal  line  exists  over  the  whole  frequency  range  of  operation. 

In  order  to  provide  shorter  interconnects  between  chips,  modern  multilayer  inte¬ 
grated  circuit  packages  for  high-performance  mainframe  computers  employ  not  only  con¬ 
ventional  striplines  but  vertical  transmission  lines  or  so-called  vias  as  well.  Because  vias 
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-may  be  of  comparable  length  with  the  striplines,  the  study  of  transient  electromagnetic 
wave  propagation  on  the  former  is  equally  indispensible  to  the  understanding  of  how  fast 
the  multilayer  integrated  circuits  can  operate.  — ^  j->  ) 

'  '  / 

Idealized  circuit  packages  consisting  of  circular  vias  going  through  circular  holes  in 
the  ground  planes  that  separate  different  layers  of  equal  thickness  are  considered.  When 
only  one  via  is  present,  we  show  that  by  treating  the  layer  between  two  ground  planes  as  a 
radial  waveguide  driven  by  a  coaxial  feed  and  making  use  of  symmetry  the  equivalent  net¬ 
work  parameters  for  each  layer  can  be  obtained.  The  multiple  layers  constitute  a  periodic 
structure  and  presents  little  difficulty  in  computing  the  total  transmission  matrix.  We  also 
study  an  infinite  array  of  vias  and  holes  with  equal  spacing  along  two  perpendicular  direc¬ 
tions.  By  virtue  of  the  image  theory,  the  problem  is  reduced  to  a  coaxial  transmission  line 
with  square  outer  wall  and  periodic  circular  diaphragm  discontinuities.  Both  variational 
approaches  as  well  as  the  integral  equation  method  can  be  applied  to  the  calculation  of  the 
characteristic  impedance,  the  junction  impedance,  and  the  effective  propagation  constant. 
The  transient  response  is  then  obtained  by  performing  the  fast  Fourier  transform  on  the 
frequency  domain  response  or  by  direct  Laplace  inversion  when  the  dispersion  of  junction 
capacitances  is  negligible.  Our  numerical  results  based  on  the  dimensions  of  an  actual 
package  show  that  the  degree  of  reflection  and  signal  distortion  becomes  intolerable  when 
the  input  rise  time  is  shorter  than  500  ps. 

The  transient  fields  of  a  current  source  on  a  layered  medium  are  calculated  using  the 
double  deformation  technique,  in  which  complex  integrals  are  deformed  in  the  transverse 
wavenumber  and  frequency  planes.  Singularities  from  thsse  complex  planes  correspond  r_ 
to  physical  modes  of  the  structure,  such  as  guided  and  leaky  waves,  and  the  relative 
importance  of  each  to  the  overall  response  can  be  discovered.  Unlike  the  Cagniard-deHoop  - 
method,  double  deformation  can  be  applied  to  dispersive  and  dissipative  media.  Also,  the  _ 


causality  of  the  electromagnetic  signal  can  be  shown  analytically. 
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A  modification  to  the  double  deformation  technique  is  presented,  which  consists  of 
splitting  the  Fourier  transform  of  the  source  current  into  two  halves,  one  for  times  before 
the  arrival  at  the  observation  point,  and  one  after.  This  greatly  increases  the  range  of 
sources  to  which  the  double  deformation  technique  can  be  applied.  Another  advantage 
of  the  modification,  the  individual  causality  and  continuity  of  each  mode,  will  be  shown. 
The  relative  importance  of  the  different  wave  modes  will  be  demonstrated,  as  well  as  the 
improvement  of  this  method  over  more  brute-force  approaches. 

Results  will  be  shown  both  for  line  and  strip  currents  on  the  surface  of  a  coated 
perfect  conductor,  for  cases  where  the  dielectric  coating  is  both  lossless  and  dissipative.  In 
most  cases,  only  a  small  number  of  modes  suffices  to  reproduce  the  important  features  of 
the  response,  including  the  arrivals  of  reflected  and  lateral  rays.  The  importance  of  each 
type  of  arrival  depends  on  certain  features  of  the  time  function,  especially  the  initial  slope. 
The  response  due  to  a  strip  current  resembles  that  of  a  line  current,  with  some  smoothing 
of  the  sharper  features. 

An  electroquasistatic  approach  in  the  spectral  domain  is  derived  to  analyze  a  pad- 
electrode  antenna  used  in  the  applications  of  low  frequency  geophysical  probing.  An  inte¬ 
gral  equation  is  derived,  which  is  then  solved  by  the  method  of  moments.  A  set  of  unit  pulse 
functions  are  used  to  represent  the  outflowing  current  distribution  on  the  pad-electrode 
surface.  The  potential  distribution  in  the  layered  medium  can  be  calculated  in  terms  of 
the  outflowing  current  distribution. 

Numerical  results  for  a  symmetrical  three-layer  medium  is  presented.  The  effects  of 
layer  thickness  and  conductivity  is  studied.  A  more  practical  measurement  environment 
with  four  different  kinds  of  conductivity  profiles  are  then  investigated.  The  effects  of  stand¬ 
off  thickness,  background  conductivity,  and  profile  depth  are  studied.  Numerical  results 
reveal  that  the  total  electrode  ament  is  sensitive  to  the  conductivity  of  the  surrounding 
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medium,  but  is  insensitive  to  the  standoff  thickness.  Thus,  it  can  be  used  to  retrieve  the 
background  conductivity  for  different  types  of  conductivity  profiles,  which  is  very  useful 
in  geophysical  probing. 

A  hybrid  method  for  radar  cross  section  calculation  of  perfect  conductors  is  dis¬ 
cussed.  The  hybrid  method  utilizes  the  electric  field  integral  equation  formulation;  it 
combines  the  method  of  moments  and  the  physical  theory  of  diffractions.  The  method  of 
moments  is  employed  to  obtain  the  matrix  equation  for  the  problems.  The  physical  the¬ 
ory  of  diffractions  is  employed  to  obtain  guess  solutions.  The  conjugate  gradient  iterative 
method  is  applied,  with  5  %  error  criterion  and  the  guess  solutions,  to  obtain  accurate 
solutions.  Numerical  results  indicate  that  the  5  %  error  criterion  and  the  guess  solution, 
provided  by  the  physical  theory  of  diffractions,  are  effective  for  obtaining  accurate  solu¬ 
tions  and  reducing  the  computation  time.  For  the  monostatic  radar  cross  section  problems, 
the  physical  theory  of  diffractions  is  used  to  provide  the  guess  solution  for  the  first  angle, 
and  for  all  subsequent  angles  the  previous  solutions  and  the  physical  optics  phase  cor¬ 
rections  are  used  to  generate  guess  solutions.  Based  on  case  studies,  it  is  found  that  the 
hybrid  method  converges  and  the  average  number  of  iterations  is  insensitive  to  the  number 
of  unknowns  and  the  shape  of  the  targets;  it  is  below  10  in  most  cases.  Therefore,  the 
computation  time  is  proportional  to  the  number  of  unknowns  to  the  second  powers. 

Cylindrical  microstrip  antennas  find  many  applications  pertaining  to  high  speed 
aircrafts  and  spacecrafts,  because  of  their  conformity  with  the  aerodynamical  structure 
of  such  vehicles.  Recently  there  has  been  some  progress  in  the  theoretical  study  of  such 
antennas,  where  the  radiation  from  various  cylindrical  microstrip  elements  was  computed 
by  assuming  an  electric  surface  current  distribution  on  the  microstrip  patch.  The  excitation 
problem  of  realizing  such  a  current  distribution  still  need  to  be  addressed.  Furthermore, 
the  input  impedance  for  the  cylindrical  microstrip  antennas  has  not  been  reported. 


We  address  the  more  realistic  problem  of  the  radiation  from  a  cylindrical  microstrip 
antenna  excited  by  a  probe.  Both  the  cylindrical-rectangular  and  the  wraparound  elements 
are  discussed.  The  current  distribution  on  the  patch  is  rigorously  formulated  using  a 
cylindrically  stratified  medium  approach.  A  set  of  vector  integral  equations  are  derived 
which  governs  the  current  distribution  on  the  patch.  This  set  of  equations  is  then  solved 
using  Galerkin’s  method  in  which  the  patch  current  is  expanded  in  terms  of  a  complete 
set  of  basis  functions  that  can  take  into  account  the  edge  singularity  condition.  The  input 
impedance  together  with  the  radiation  pattern  are  derived  both  exactly  and  in  the  small 
substrate  thickness  limit  where  a  single  mode  approximation  is  employed. 

A  full-wave  analysis  is  developed  to  rigorously  study  the  patch  current  distribu¬ 
tion,  the  radiation  fields,  and  the  input  impedance  for  the  wraparound  and  rectangular 
microstrip  antennas  excited  by  a  current  probe.  The  analysis  is  valid  for  both  thick  and 
thin  substrates. 

For  thick  substrates,  hybrid  modes  are  excited,  and  only  in  the  case  of  axially  sym¬ 
metric  modes  (»  =  0)  the  TEom  is  decoupled  from  the  TMom»  and  modes  of  different 
parity  do  not  couple.  In  the  wraparound  patch,  all  current  modes,  with  no  axial  varia¬ 
tion,  radiate  weakly,  and  consequently,  they  have  a  narrow  bandwidth.  The  presence  of 
the  dielectric  substrate  widens  the  bandwidth  and  broadens  the  radiation  pattern.  The 
radiation  pattern  is  insensitive  to  the  substrate  thickness  (especially  for  high  dielectrics). 
The  rectangular  patch  is,  generally,  slightly  less  radiating  than  the  wraparound. 

A  rigorous  analysis  of  the  resonant  frequency  problem  of  both  the  cylindrical-  rect¬ 
angular  and  the  wraparound  microstrip  structures  is  presented.  The  problem  is  formulated 
in  terms  of  a  set  of  vector  integral  equations.  Using  Galerkin’s  method  is  solving  the  in¬ 
tegral  equations,  the  complex  frequencies  are  studied  with  sinusoidal  basis  functions.  The 
effect  of  the  edge  singularity  on  the  convergence  is  invesgated.  A  perturbative  approach  is 
also  used  to  compute  the  complex  resonant  frequencies  in  the  thin  substrate  limit. 


In  microelectronic  computer  packaging,  a  problem  of  practical  interest  is  the  study 
of  propagation  characteristics  of  a  shielded  microstrip  line  in  the  presence  of  crossing  strips. 

The  microstrip  line  and  the  crossing  strips  are  placed  at  two  different  interfaces  of 
a  three-layer  medium  bounded  by  two  parallel  conducting  plates  as  shown  in  Fig.  1.  The 
crossing  strips  are  assumed  to  be  periodic  and  the  three  layers  to  be  uniaxially  anisotropic. 

The  network  analytical  method  of  electromagnetic  fields  is  applied  for  the  hybrid¬ 
mode  analysis  of  the  frequency  dependent  characteristics  of  the  structure.  The  transverse 
fields  in  each  region  are  expressed  by  their  Fourier  transform  in  the  x-direction  and  a 
Floquet  harmonic  representation  in  the  y- direction.  Using  Maxwell’s  equations  and  ap¬ 
plying  boundary  conditions,  we  obtain  a  set  of  coupled  integral  equations  for  the  current 
distributions  /,,  and  J„  on  the  microstrip  line  and  the  crossing  metallic  strips,  respectively. 

The  determinantal  equation  for  the  dispersion  relation  can  be  derived  by  applying 
Galerkin’s  procedure  to  the  derived  set  of  coupled  integral  equations.  The  Brillouin  dia¬ 
grams  are  obtained  numerically  by  seeking  the  roots  of  the  obtained  eigenvalue  equation. 
The  stop-band  properties  are  numerically  presented  as  a  function  of  the  spacing,  length, 
and  width  of  the  crossing  metallic  strips.  The  effects  of  the  material  and  geometrical 
parameters  are  also  investigated. 

We  analyse  the  quasi-electrostatic  fields  generated  by  an  electrode  mounted  on  a 
perfectly  conducting  pad  of  a  finite  extent  in  a  planar  stratified  medium.  The  axis  of 
stratification  is  perpendicular  to  the  pad.  The  electrode  Is  maintained  at  a  prescribed 
potential  relative  to  the  pad.  This  problem  serves  as  a  canonical  problem  which  is  useful 
to  model  some  of  the  antennas  employed  in  geophysical  applications.  It  is  more  pertinent 
to  tools  which  are  mounted  on  a  pad  and  operating  at  low  frequencies. 
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We  derive  an  integral  equation  in  the  spectral  domain  for  the  normal  current  dis¬ 
tribution  on  the  pad’s  surface.  A  method  of  moments  is  then  applied  to  solve  the  integral 
equation.  A  proper  set  of  local  basis  functions  is  used  to  represent  the  unknown  normal 
current  distribution  on  the  pad.  The  basis  functions  are  chosen  to  satisfy  the  edge  con¬ 
dition  as  accurately  as  possible  so  that  the  current  distribution  can  be  represented  by  a 
reasonable  number  of  basis  functions.  Both  the  current  distribution  on  the  pad’s  surface 
and  the  current  flow  pattern  in  the  stratified  medium  are  computed. 

The  exact  image  method,  previously  introduced  for  problems  involving  sources  in 
real  space  with  two  homogeneous  media  and  a  planar  interface,  is  extended  to  problems 
involving  sources  in  complex  space.  The  idea  seems  straightforward  because  the  analytic 
functions  of  real  space  source  position  can  be  analytically  extended  to  functions  of  complex 
space  source  position,  which  has  been  applied  before  with  success  in  diffraction  problems. 
The  complex  source  point  theory  gives  a  possibility  to  analyse  Gaussian  beam  problems 
with  a  simple  source  instead  of  a  complicated  source  in  real  space.  This  has  been  tested 
by  analysing  asymptotic  (far  field)  expressions  for  a  reflecting  and  a  transmitted  Gaussian 
beam  resulting  in  well  known  Goos-Hanchen  and  angular  shifts  for  the  reflected  beam  and 
apparent  position  for  the  transmitted  beam.  The  present  method  suggests  itself  for  more 
exact  calculation  of  the  fields  in  these  beams. 

The  purpose  of  this  paper  is  to  analyse  the  transient  behavior  of  nonuniformly 
coupled  and/or  frequency  dependent  transmission  line  systems  terminated  with  nonlinear 
loads.  In  the  case  of  frequency  independent  line  systems,  the  transient  analysis  is  directly 
performed  in  the  time  domain  based  on  variable  transformations  and  the  method  of  charac¬ 
teristics.  When  the  system  is  frequency  dependent,  a  new  hybrid  method  of  combining  the 
frequency  and  time  domain  analysis  is  presented  for  dealing  with  the  nonlinear  problem. 

Theory  for  quasi- TEM  modes  propagating  in  a  transversely  inhomogeneous  (  multi¬ 
dielectric)  longitudinally  uniform  transmission  line,  is  derived  for  transient  signals.  It  is 


seen  that,  while  the  starting  point  for  the  theory  is  completely  different,  the  result  is  similar 
to  the  time-harmonic  theory,  and  previously  derived  properties  for  propagating  modes  also 
apply  in  the  transient  case. 

The  analysis  of  resonance,  input  impedance,  and  radiation  of  the  elliptic  disk  mi¬ 
crostrip  structure  is  rigorously  formulated,  using  the  scalar  and  vector  Mathieu  transforms. 
With  the  help  of  these  transforms,  the  resonance  frequencies  of  the  structure  can  be  de¬ 
rived  exactly  using  Galerkin’s  method  and  approximately  using  a  perturbational  approach. 
Expressions  for  the  input  impedance  and  the  radiation  pattern  are  also  obtained. 

The  transient  electromagnetic  radiation  by  a  vertical  electric  dipole  on  a  two-  layer 
medium  is  analysed  using  the  double  deformation  technique,  which  is  a  modal  technique 
based  on  identification  of  singularities  in  the  complex  frequency  and  wavenumber  planes. 
Previous  application  of  the  double  deformation  technique  to  the  solution  of  this  problem 
is  incomplete  in  the  early  time  response.  We  show  that  the  existence  of  a  pole  locus  on 
the  negative  imaginary  frequency  axis,  which  dominates  the  early  time  response,  proves 
crucial  in  obtaining  the  solution  for  all  times.  A  variety  of  combinations  response,  proves 
crucial  in  obtaining  the  solution  for  all  times.  A  variety  of  combinations  of  parameters 
are  used  to  illustrate  the  double  deformation  technique,  and  results  will  be  compared  with 
those  obtained  via  explicit  inversion,  and  a  single  deformation  method. 

Many  integrated  circuits  contain  strip  lines  at  different  heights  that  run  parallel  or 
perpendicular  to  each  other.  And  we  have  investigated  reliable  models  for  these  struc¬ 
tures.  First  the  capacitances  associated  with  two  offset  parallel  strips  at  different  heights 
between  ground  planes  are  computed  using  the  conformal  mapping  approach.  As  an  ex¬ 
tension,  a  simplified  circuit  of  parallel-plate  lines  with  transverse  ridges  is  introduced  to 
model  two  parallel  strips  with  perpendicularly  crossing  strips  on  top.  We  treated  it  as 
a  distributed  circuit  consisting  of  transmission  lines  segments  with  periodical  capacitive 
loading.  In  order  to  calculate  the  coupling  between  two  lines,  we  reduced  this  structure 
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to  two  equivalent  single  line  circuits,  viz.  the  even  mode  and  the  odd  mode  circuits.  The 
Laplace  transform  approach  can  be  easily  applied  to  find  out  the  transient  response.  The 
numerical  computation  carried  out  for  various  environments  shows  that  the  crossing  strips 
will  cause  serious  trouble  for  signals  with  a  rise  time  of  less  than  50ps  to  propagate  along 
a  distances  of  2cm  or  longer. 

Basically,  vias  in  a  multilayered  integrated  circuits  are  treated  like  transmission 
lines  with  loadings  where  they  encounter  holes  in  ground  planes  separating  different  lay¬ 
ers.  We  have  modeled  a  ground  plane  with  a  hole  and  a  circular  conductor  at  the  center 
of  the  hole  as  a  radial  waveguide,  which  in  turn  is  connected  to  the  via,  another  section 
of  transmission  line.  Thus  by  computing  the  characteristic  impedance  of  the  former,  we 
have  derived  the  equivalent  load  impedance  of  the  via  hole.  The  load  impedance  is  one 
important  parameter  in  calculating  the  transient  propagation  along  vias. 

The  scattering  of  electromagnetic  waves  from  a  randomly  perturbed  peridic  surface 
is  solved  using  the  Extended  Boundary  Condition  (EBC)  method.  The  scattering  from 
periodic  surface  is  solved  exactly  using  the  EBC  method  and  this  solution  is  used  in  the 
small  perturbation  method  to  solve  for  the  scattered  field  from  a  randomly  perturbed  peri¬ 
odic  surface.  The  random  perturbation  is  modeled  as  a  Gaussian  random  process  and  the 
surface  currents  and  the  scattered  fields  are  expanded  and  solved  up  to  the  second  order. 
The  theoretical  results  are  illustrated  by  calculating  the  bistatic  and  backscattering  coef¬ 
ficients.  It  is  shown  that  as  the  correlation  length  of  the  random  roughness  increases,  the 
bistatic  scattering  pattern  of  the  scattered  fields  show  several  beams  associated  with  each 
Bragg  diffraction  direction  of  the  periodic  surface.  When  the  correlation  length  becomes 
smaller,  then  the  shape  of  the  beams  become  broader.  The  results  obtained  using  the  EBC 
method  is  also  compared  with  the  results  obtained  using  the  Kirchhoff  approximation.  It 
is  shown  that  the  Kirchhoff  approximation  results  show  quite  a  good  agreement  with  EBC 
method  results  for  the  VV  and  HH  polarized  backscattering  coefficients  for  small  angles 


of  incidences.  However,  the  Kirchhoff  approximation  does  not  give  depolarized  returns 
in  the  backseat tering  direction  whereas  the  results  obtained  using  the  EBC  method  give 
significiant  depolarized  returns  when  the  incident  direction  is  not  perpendicular  to  the  row 
direction  of  the  periodic  surface. 

The  analysis  of  resonance,  input  impedance  and  radiation  of  the  elliptic  disk,  mi- 
crostrip  structure  is  rigorously  formulated,  using  the  Scalar  and  Vector  Mathieu  Trans¬ 
forms.  With  the  help  of  these  transforms,  the  resonance  frequencies  of  the  structure  can 
be  derived  exactly  using  Galerkin’s  method  and  approximately  using  a  perturbational  ap¬ 
proach.  Expressions  for  the  input  impedance  and  the  radiation  pattern  are  also  obtained. 

Theory  for  quasi-TEM  modes  propagating  in  a  transversely  inhomogeneous  (multi¬ 
dielectric)  longitudinally  uniform  transmission  line,  previously  derived  for  time-harmonic 
waves,  is  derived  for  transient  signals.  It  is  seen  that,  while  the  starting  point  for  the  theory 
is  completely  different,  the  result  is  similar  to  the  time-harmonic  theory,  and  previously 
derived  properties  for  propagating  modes  also  apply  in  the  transient  case.  The  range  of 
applicability  is  discussed  with  a  simple  example. 

Exact  image  method,  recently  introduced  for  the  solution  of  electromagnetic  field 
problems  involving  sources  above  a  planar  interface  between  two  homogeneous  media,  was 
originally  applied  to  a  restricted  pair  of  medium  parameters  to  obtain  a  well  behaved  image 
located  in  the  ’proper’  complex  half  space.  It  is  demonstrated  here  with  an  example  that 
the  ’proper’  half  space  limitation  can  be  released  to  increase  the  applicability  of  the  exact 
image  theory.  However,  it  is  shown  that  for  certain  media,  numerical  difficulties  in  the 
field  integrals  may  be  encountered,  due  to  crossing  of  branch  cut  lines  on  the  complex 
integration  plane.  This  may  occur  when  the  medium  where  the  field  is  to  be  calculated 
is  more  lossy  than  the  other  medium.  Methods  for  the  image  integration  to  obtain  better 
convergence  for  more  general  media,  are  discussed. 
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Simple  approximation  for  diffraction  surface  currents  on  a  conducting  half  plane,  due 
to  an  incoming  plane  wave,  is  given  in  terms  of  a  line  current  (monofile)  in  complex  space. 
When  compared  to  the  approximation  by  a  current  located  at  the  edge,  the  diffraction 
pattern  is  seen  to  improve  by  an  order  of  magnitude  for  a  minimal  increase  in  computation 
effort.  Thus,  the  inconvenient  Fresnel  integral  functions  can  be  avoided  in  quick  calcula¬ 
tions  of  diffracted  fields  and  the  accuracy  is  seen  to  be  good  in  other  directions  than  along 
the  half  plane.  The  method  can  be  generally  applied  to  problems  involving  planar  metal 
edges. 

We  have  also  proposed  a  way  to  account  for  the  effect  of  complicated  geometry  from 
the  point  of  view  of  continuous  transmission  line  model  by  letting  the  coupling  between 
parallel  lines  in  multilayered  integrated  circuits  to  be  nonuniform.  Previously,  we  devised 
a  scheme  that  combines  the  method  of  characteristics  and  perturbational  series  to  simplify 
the  computation  of  the  transient  response  from  the  coupled  transmission  line  equations.  A 
new  transformation  for  decoupling  now  enables  us  to  generalize  this  formulation  to  calcu¬ 
late  the  near-end  and  far-end  crosstalks  to  very  high  accuracy,  given  arbitrary  positional 
dependence  for  both  capacitive  and  inductive  coupling  coefficients. 

The  transient  response  of  fundamental  sources,  such  as  dipole  and  line  current,  was 
carefully  analyzed.  With  the  double- deformation  technique,  which  is  a  modal  technique 
based  on  identification  of  singularities  in  the  complex  frequency  and  wave  number  planes, 
we  are  able  to  obtain  both  early  and  late  time  response  very  efficiently.  Some  results  for 
vertical  electric  dipole  excitation  on  a  two-layer  medium  have  been  published.  Recently, 
we  have  discovered  a  general  scheme  of  breaking  up  the  integrands  so  that  sources  with 
arbitrary  time  and  space  dependence  can  be  incorporated  into  our  formulation  without 
sacrificing  convergence. 
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The  exact-image  method,  recently  introduced  for  the  solution  of  electromagnetic  field  problems  involving  sources 
above  a  planar  interface  between  two  homogeneous  media,  is  shown  to  be  valid  also  for  sources  located  in  complex 
space,  which  makes  its  application  possible  for  Gaussian-beam  analysis.  It  is  demonstrated  that  the  Goos- 
Hinchen  shift  and  the  angular  shift  of  a  TE-po lamed  beam  are  correctly  given  as  asymptotic  reeulta  by  the  exact- 
reflection-image  theory.  Also,  the  apparent -image  location  giving  the  correct  Gaussian  beam  transmitted  through 
the  interface  is  obtained  as  another  asymptotic  check.  The  theory  described  here  makee  it  possible  to  calculate  the 
exact  coupling  from  the  Gaussian  beam  to  the  reflected  and  refracted  beams  as  well  as  to  the  surface  wave. 


INTRODUCTION 

The  Gaussian  beam  is  an  important  model  for  radiation 
from  large-aperture  antennas  and  laser  sources,  from  which 
the  energy  is  sent  in  a  relatively  narrow  space  angle.  The 
introduction  by  Deschamps1  of  a  complex-space  dipole  to 
represent  mathematically  a  source  for  the  Gaussian  beam 
made  it  possible  to  extend  analytical  results,  derived  earlier 
for  real-space  sources  to  Gaussian-beam  excitations.  The 
basic  Gaussian  beam  originates  from  a  dipole  in  complex 
space,  whereas  multipoles  of  higher  order  can  be  shown  to 
produce  more-complex  Hermite-Gaussian  and  Laguerre- 
Gaussian  beams.2'4 

The  basic  problem  in  Gaussian-beam  analysis  of  a  planar 
interface  of  two  media  is  that  of  reflection  and  refraction  of 
the  beam  at  the  interface.  This  problem  has  been  treated  by 
using  a  complex -space  dipole  method8  and  by  other  meth¬ 
ods.8-8  The  general  field  expressions  in  the  exact  formula¬ 
tions  are  involved,  but  asymptotic  considerations  for  narrow 
beams  permit  a  simpler  analysis  of  reflection  and  transmis¬ 
sion  problems  and  lead  to  the  well-known  shift  phenomena 
of  the  reflecting  beam.  The  most  famous  of  these  was  first 
demonstrated  by  Goos  and  H&nchen®  in  1947  and  analyzed 
by  Artman10  in  1948  and  involved  a  parallel  shift  of  the  beam 
reflected  from  an  interface  in  the  denser  of  two  lossless 
dielectric  media.  Other  shifts  that  can  be  defined  for  nar¬ 
row  beams  are  the  angular  shift8-11  and  the  focal  shift.12 
Infinities  arising  in  the  idealized  theory  for  the  Goos-Hin- 
chen  shift  can  be  avoided  by  using  a  more-realistic  analy- 
sis.8t3,u 

The  purpose  of  the  present  paper  is  to  give  a  more-general 
account  of  the  exact-image  theory,  which  was  introduced 
elsewhere  for  sources  in  the  air  above  a  planar  interface  of  a 
homogeneous  medium1818  (the  Sommerfeld  problem),  and 
to  allow  the  original  source  to  be  in  complex  space.  Thus  the 
theory  is  applicable  to  global  Gaussian-beam  analysis  in 
exact  form.  In  comparison  with  other  exact  methods  based 
on  Sommerfeld  integrals  (inverse-Fourier-transform  inte¬ 
grals  for  the  fields),  the  present  theory,  which  deals  with 
equivalent  sources  and  their  integration,  gives  more  physical 
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insight  and  does  not  resort  to  special  integration  procedures 
that  are  dependent  on  the  field  point  For  the  application  of 
the  exact-image  theory,  functions  characterizing  the  image 
sources  are  needed,  and  methods  for  their  computation  are 
presented  in  Refs.  15  and  16.  These  functions  need  to  be 
calculated  only  once  in  the  computer  memory,  after  which 
the  field  calculation  involves  converging  well-behaved  inte¬ 
grals.  The  exact  image  of  a  point  source  is  a  line  source, 
which  in  general  lies  in  complex  space.  To  check  the  theory, 
known  asymptotic  expressions  for  the  reflected  beam  shifts 
and  the  apparent-transmission-image  location  are  shown  to 
arise  as  special  cases  when  the  exact  line  image  is  approxi¬ 
mated  by  a  point  image. 


REFLECTION-IMAGE  THEORY 

The  notation  applied  here  is  based  on  that  used  in  Ref.  15 
except  that  we  consider  two  half-spaces  with  parameters 
mimo  and  1 1(0  for  &  •  f  >  0  and  parameters  M2M0  and  <2<o  for  G  •  t 
<  0.  In  Refs.  15  and  16  the  medium  fi  •  f  >  0  was  assumed  to 
be  air,  but  generalization  for  any  medium  pair  is  straightfor¬ 
ward  if  instead  of  k  we  write  «  hG»i«i)1/2,  where  k  m 
<*>(ao«o)1/2  and  we  understand  that  a  “  mj/mi  and  tj /ej  every¬ 
where.  To  avoid  unnecessary,  although  not  excessive,  com¬ 
plication,  the  theory  here  is  limited  to  dielectric  media  only 
with  at  »  mj  ■  1. 

As  another  generalization,  the  original  source  is  taken  to 
be  in  complex  space.  Because  the  theory  in  Refs.  15  and  16 
was  analytically  derived  and  no  use  of  the  tacit  assumption 
of  a  real  source  location  was  made,  the  final  expressions  can 
also  be  applied  for  complex  source  locations.  The  same  idea 
was  applied  earlier  with  success  when  field  expressions  that 
were  derived  for  problems  with  sources  in  real  space  were 
generalized  for  sources  in  complex  space,17  resulting  in  solu¬ 
tions  for  problems  with  Gaussian-beam  excitation. 

Let  us  consider,  for  simplicity,  the  source  «?(f)  *  0lL6{t  - 
iCih  +  >6)J,  which  is  a  dipole  with  the  direction  of  the  unit 
vector  0  and  the  location  at  the  real  height  h  from  the  inter¬ 
face  d  • )  ■  0  and  the  imaginary  depth  j6.  In  the  calculation 
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of  the  fields  in  the  half-space  1,  the  half-space  2  can  be 
replaced  by  the  exact-image  source,  as  was  shown  in  Ref.  15: 

*7,0.  P)  -  ~  [/,(P>  +  *+<P)]  M  •  Je(r) 

~a-^—f,(p)a-v[vt.jcm].  d) 

"i  * 


Here,  p  is  an  integration  variable,  and  the  image  function 
f,(p)  can  most  conveniently  be  calculated  through  the  fol 
lowing  Bessel-function  series16: 


/,(P) 


— 8<  V  /«  —  1\" 

«2  —  1“J\<  +  1/  P  ' 


fliP)  * 


-2J2(p) 

P 


(2) 


Further,  4+(p)  denotes  6(p  -  0+),  t  denotes  a  component 
transverse  to  Ci,  and  c  denotes  the  reflection  operation 


d,  -  C  •  d,  JC(M  «  C  •  <?(C  •  P),  C  -  I  -  2uu.  (3) 


The  field  from  the  image  source  can  be  written  as  a  fourfold 
integral  over  space  and  the  parameter  p: 


and  lie  parallel  to  the  complex  z  plane,  with 

argtz  +  h  +  jd  •  6)  »  argOB).  (8) 

Requiring  that  the  additional  condition  Re[z)  *•  0  be  valid, 
the  complex  distance  function  satisfies  D  *  0  for  all  field 
points  in  Re(z)  >  0,  and  the  field  integrand  in  Eq.  (4)  is 
nonsingular.16  This  condition  defines  the  branch  of  B,  or 
the  square  root  («2  -  €,)1/J  through 

R*l/<«;  -  e:>,/2]  *  Re{(<,  -  f2)1/2)  >  0,  (9) 

implying  for  example,  that  if  n,  »2  “  eet  <  f,  are  real,  then 
U  -  l)1/J  ■  -/( 1  -  provided  that  (1  -  «)1/J  >  0.  For 
some  combination  of  lossy  medium  parameters  there  may  be 
uncertainty  in  choosing  the  correct  branch  of  the  distance 
function  D  to  obtain  the  best  convergence.18  A  close  study 
of  the  integration  path  on  the  complex  z'  plane  will  reveal 
that  a  branch-cut  line,  starting  from  the  branch  points  at  z' 
*  z  ±  j(>,  defining  the  converging  branch  of  the  Green  func¬ 
tion,  may  be  crossed  for  some  combination  of  parameter 
values  (]  and  «2  when  p  moves  from  0  to  ®.  In  this  case  it  is 
possible  to  take  the  image-current  line  in  the  wrong  half¬ 
space  with  Re[z']  >  0,  which  makes  the  integrand  converging 
again.  This  question  is  the  subject  of  a  forthcoming  paper. 


oj  |  Gi(D)  •  p)d V'dp,  (4) 


with 


G,(D) 


exp(-jk,Dl 
4  tD 


D(K  r.  p)  »  [(?  -  r  +  iip/jB )  -  v  -  r  +  apljB) ]1/2, 

B  -  kU2  -  «,)‘*  (5) 


The  image-current  expression  contains  Bessel  functions 
</„(p),  which  are  convergent  only  if  the  integration  parame¬ 
ter  p  is  real.  For  a  point  source,  Eq.  ( 1 )  can  also  be  written  so 
that  the  integration  parameter  p  is  absent,  because  in  the 
field  integral  [Eq.  (4)),  the  coordinate  z  ■  d  •  t  and  p  are 
related  through  the  distance  function  D  in  the  argument  of 
the  Green  function.  The  expression  [Eq.  (1)]  can  thus  be 
written 

Jt{r)  *  -1L  *  — -j  u(u  •  0)4,.(f  +  uh  —  j6c)  -  jk^t  -  1  )inIL 

X  {^.IpUMp  ~  jS,)  -  m  ■  0)f,[p(t)]i(p  -  j6t) 

-  Cl  —  f\p(z)\&  vm  +  ah-  >6r>| . 

(6) 


with 


pit)  ■  -jB(z  +  h  +  ja-  6).  (7) 

To  obtain  an  exponentially  converging  Green  function,  we 
must  select  the  branch  of  the  distance  function  D  so  that 
Im(k|Dj  <  0.  Also,  to  obtain  a  converging  image-current 
function,  the  path  of  z  integration  must  be  chosen  so  that  p 
in  Eq.  (7)  is  real  and  positive,  in  which  case  all  the  Bessel 
functions  in  Eqs.  (2)  converge.  This  means  that  the  image 
line  must  start  at  the  point  A  *  -uh  +  Sc  in  complex  space 


REFLECTED  GAUSSIAN  BEAM 

A  unique  global  definition  for  the  Gaussian-beam  field  does 
not  seem  to  exist;  instead,  the  term  Gaussian  beam  is  under¬ 
stood  as  an  asymptotic  property  of  radiation  fields  close  to 
the  axis  of  the  main  radiation.  Thus  many  fields  with  the 
same  asymptotic  quadratic  exponential  behavior  are  called 
Gaussian  beams.  (The  present  method  is  not,  however, 
dependent  on  the  definition  of  the  Gaussian  beam,  since  it 
works  for  any  original  sources.)  One  example  is  the  field 
from  a  point  source  in  complex  space  with  the  imaginary 
part  of  the  position  vector  6  *  Ub  cos  8  -  Uyb  sin  8. 

The  radiation  beam  of  the  point  source  is  obtained  in  the 
direction  of  -S.  The  reflected  beam  is  obtained  from  the 
image  source  [Eq.  (6)],  and.  for  points  far  enough  from  the 
interface,  the  field  can  be  approximated  by  a  Gaussian  beam 
if  the  image  source  can  be  approximated  by  a  point  source. 
Of  course,  the  correct  field  is  obtained  also  in  space  points, 
where  it  cannot  be  described  as  a  Gaussian  beam.  To  dem¬ 
onstrate  the  validity  of  the  exact-image  theory  for  this  prob¬ 
lem.  we  show  that  the  Goos-Hanchen  shift  and  the  angular 
shift  of  the  reflected  beam  are  obtained  as  asymptotic  re¬ 
sults  from  the  exact-image  field  expression.  To  keep  the 
notation  simple,  let  us  consider  a  two-dimensional  problem 
with  the  line  source  J(f)  *  uji(y  +  jb  sin  8)4(z  -  h  —  jb 
cos  8),  where  8  is  the  angle  between  a  and  6.  This  assump¬ 
tion  simplifies  the  image  source  [Eq.  (1)]  into 

p)  -  /,(p)Jr(A)  -  +  jb  sin  8) 

X  Hz  -t-  h  +jb  cos  8),  (10) 

where  the  function  f\(p)  is  that  given  in  Eqs.  (2).  The  Green 
function  [Eqs.  (5)]  integrated  in  the  x  direction  produces  the 
two-dimensional  Green  function  in  the  medium  1, 

G,(D)-^I  +  ^-2rvj  G,(D),  GX(D)  -  J-  W0I2,(*,D),  <H) 
where  H0l2’l*l  *®  the  Hankel  function. 
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The  field  calculated  from  the  exact-image  source  is  valid 
in  any  point  in  the  upper  half-space,  but  since  we  are  inter¬ 
ested  in  Gauss ian-beam  properties,  let  us  consider  the  far 
field  with  lAtDl  »  1,  in  which  case  the  image  can  be  approxi¬ 
mated  by  a  point  source.  Also,  the  asymptotic  expression 
for  the  Hankel  function  can  now  be  used: 

H0,2](k{D)  *  ^  exp <-;*,D).  (12) 

D  here  represents  the  complex  distance  between  the  field 
point  f  and  the  complex  integration  point  f'  »  jup! 
fci  (<  -  l)l/l,  where  P  *  -ud  +  j6c,  if  p  integration  is  done 
separately.  Because  the  image  function  flip)  is  decaying, 
the  effective  p-integration  range  can  be  regarded  as  small 
with  respect  to  If  -  H  and  we  can  write 


D  »  O0  -  D0  -  [<f  -  n  ■  if  -  n)u\ 

*1 


u- Lb 

(r - 1)1/2  ’ 


w 


Dn 


(13) 


By  applying  expression  (12),  we  earn  write  the  field  integral 
on  the  plane  x  *  0  as 

£{»  *  -JwhquJG^Dq)  I"  /.(ple'^dp.  (14) 

Jo 

If  expression  (14)  can  be  written  in  the  form 


with  the  branch  of  the  square  root  so  defined  that 
Im[(*  -  sin2  0)l;'2]  <  0.15  Equation  (19)  explains  the  Goos- 
Hanchen  shift  and  the  angular  shift  when  the  Taylor-ap- 
proximation  condition  is  valid.  In  fact,  because  w0  ■  -6c/fe, 
we  have  q0~  cos  0/(«  -  l)l/2,  and  Eq.  (17)  represents  the  TE 
reflection  coefficient  of  a  plane  wave  coming  at  the  angle  0. 
Let  us  consider  the  two  possibilities  with  real  r: 


(1)  a  <  sin2  9,  which  can  only  happen  for  <,  >  «2.  In  this 
case,  3  is  real,  and  thus  the  real  pan  of  the  location  vector  is 
shifted  by 


3 


2  u 

fe^sin2  9  -  t)l/J 


(20) 


which  means  a  parallel  shift  of  the  reflecting  beam  (Fig.  1). 
This  is  the  well-known  Goos-Hanchen  shift.  Because  s  is  in 
the  -u  direction,  the  parallel  shift  is  s  sin  9. 

(2)  (  >  sin2  9.  In  this  case  the  square  root  is  real  and  the 
shift  3  is  imaginary  (I  ■  jsu),  which  means  that  the  image 
line  is  shifted  an  imaginary  distance  in  the  &  direction  from 
the  original  location.  Because  the  imaginary  part  of  the 
position  vector  determines  the  direction  of  the  beam,  the 
beam  is  shifted  angularly  from  its  mirror-image  direction  w0 
to  the  direction  determined  by  6C  -  j3  (Fig.  2).  If  the  shift 
angle  A 9  is  small,  it  satisfies 


tan(A0)  *  —  cos  9  « 
o 


2  cos  9 

k,6(f  -  sin20)I/s 


(21) 


£(3)  *  -  ;W0«,fGi(£,o)  expi  jkxu  ■  3)  expi-jkxw0  ■  s ) 

X  j  fxip)  exp(-pq0)dp,  (15) 
Jo 


in  the  vicinity  of  the  radiating  direction  w  *  w0,  where  wo 
corresponds  to  a  field  point  f0  on  the  axis  of  the  reflected 
beam,  the  field  can  be  thought  of  as  arising  from  the  current 
line  shifted  by  the  vector  5  from  the  mirror-image  location  f. 
Approximation  (15)  is  obviously  possible  if  expression  (14) 
can  be  expanded  as  a  power  series  in  terms  of  the  small 
difference  vector  ill  -  t&o.  This  is  possible  only  for  sufficient¬ 
ly  narrow  beamwidths.  After  some  steps  of  Taylor  expan¬ 
sion,  the  expression  for  3  can  be  written  in  the  form 


a- . - 

p/,(p)  exp(-pq0)dp 

) 

A,(«  -  l)l/2 

J ^  A  (P>  exp(-p<70)dp 

(16) 


where  qom  &•  il>o/(«  ~  l)l/l-  To  obtain  an  expression  for  i, 
the  following  integral  identities  are  needed: 


W*dP-'>--(g2+--l-)1- 

1  <?  +  (<72  +  l)l/2 


(17) 


'  pU(p)t-»Ap  .  -2  1  -  <gl ♦  V-2  .  (18) 

'  P  iq1  +  1>1/2  q  +  «72  +  l)1'* 


Thus  we  have 


2 ju 


2ju 


*,(f  -  1),'2(<J02  +  1)1/2  fe[«2  -  «,  +  «,(u  •  ii0)2|l/2 

2ju 

Jt,(»  -  sin20)l/2 


The  previous  expressions  are  valid  only  for  sufficiently 
narrow  beams,  a  condition  that  presumes  sufficiently  large 
values  for  b.  They  are  not  valid  if  the  Taylor  expansion  is 
not  applicable,  which  happens  at  and  close  to  the  branch 
point  of  the  square-root  function  (q02  +  1)1/2,  i.e.,  at  sin2  9  * 
sin2  9e  *  e,  where  Eq.  (19)  would  predict  an  infinite  shift. 
This  is  the  definition  of  the  critical  angle  9C,  for  which 
(<?o2  +  1)1/2  *  0.  For  angles  0  *  0C  we  can  write,  from 
expression  (14)  and  Eq.  (17), 

£(r)  *  -jvtnMJGx(DQ)  exp(-2;[2(coe  0  -  cos0c)/(l  -  t)l/2]l/2l, 

(22) 


Fig.  1.  Goos-Hanchen  shift  of  the  Gaussian  beam  arising  from  an 
approximate-image  point  source  with  a  real  shift  vector  1  from  the 
mirror-image  point  -uh  +  j6c. 


■  < 


i 


0 


(19/ 
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Fig.  2.  Angular  shift  of  the  Gaussian  beam  arising  from  an  approxi¬ 
mate-image  point  source  with  an  imaginary  shift  vector  i  from  the 
mirror-image  point  -Cih  +  jS,. 


which  cannot  be  described  by  a  simple  shift  of  the  image 
source  because  the  reflecting  beam  is  distorted  and  not  ex¬ 
actly  of  Gaussian  form.  To  find  out  the  shift  of  the  maxi¬ 
mum  of  the  reflecting  beam,  some  numerical  analysis  must 
be  made,  as  in  Refs.  6, 13,  and  14. 

For  complex  t  values,  there  are  both  real  and  imaginary 
parts  for  the  vector  8,  .hich  means  a  combined  Goos-Han- 
chen  and  angular  shift.  In  this  case,  the  beam  does  not  enter 
at  the  critical  angle,  because  it  is  now  complex.  For  beams 
sufficiently  narrow  the  real  and  imaginary  parts  of  i  in  Eq. 
(19)  determine  the  Goos-Hanchen  and  angular  shifts,  re¬ 
spectively. 

This  simple  test  serves  as  an  asymptotic  verification  of  the 
exact-image  theory  for  the  Gaussian-beam  reflection  analy¬ 
sis.  The  field  of  the  reflected  beam  can  be  calculated  at  any 
point  from  the  exact  formulation  [Eq.  (4)|  and  also  at  points 
where  it  cannot  be  described  as  a  Gaussian  beam  and  where 
the  asymptotic  formulas  are  not  valid.  This  also  includes 
the  coupling  to  the  surface  wave,  which,  however,  is  small 
when  8  is  not  *0. 

TRANSMISSION-IMAGE  THEORY 

The  image  theory  for  fields  transmitted  through  a  planar 
interface  was  developed  in  Ref.  16,  and  the  generalization 
discussed  for  reflection-image  theory  also  applies  here.  To 
summarize  the  result,  the  upper  dielectric  half-space  is  re¬ 
placed  by  the  transmission- image  source: 

J,(F  p)  -  Ciu  ■  J(8)  Mp)  +  p) 

+  J,(f)[8+(p)  +  F{'(Bz,  p)]  +  u[F,(Bz,  p) 

-  F.lBz.plJH'U.p^.JlH,  (23) 

Here,  the  original  source  is  <7(^),  and  z  denotes  its  coordinate. 
Further,  we  have 


fl  “  *,(<  -  1)1/2  *  *(«j  -  r,)172, 


H(»,p)-[rJ  +  (p/B)J]1/s, 


F.U.p)  -  J0(p) 


(24) 

(25) 


X 


V  [ill  1  ~  [(p/r)2  +  x]1/21 m 

",  1*  +  1  1  +  l(p/r)2  +  l]l/2J 


Jtmip).  (26) 


and  the  primes  in  functions  H  and  F  denote  differentiation 
with  respect  to  p. 

The  field  is  obtained  from  an  integral  similar  to  that  in  Eq. 
(4),  but  the  medium  is  now  «j.  Also,  the  distance  function  D 
is  here  defined  to  be 


do,  r ,  P)  -  i[f  -  p'  -  amx\  p»  •  [f  -  p-  -  amz\  p)]il/2. 

(27) 

Again,  to  obtain  a  converging  image  function,  the  parame¬ 
ter  p  must  be  real.  To  obtain  a  converging  Green  function, 
Im[AjZ)|  must  be  nonpositive,  which  defines  the  branch  of 
the  D  function.  Also,  in  this  case,  for  certain  lossy  media, 
the  branch  cut  of  the  Green  function  may  be  crossed  if  the 
image  is  restricted  to  the  half-space  Re[z']  >  0,  leading  to  a 
nonconverging  Green  function.  This  problem  can  be  over¬ 
come  by  taking  the  other  branch  of  the  H(z,  p)  function. 


TRANSMITTED  GAUSSIAN  BEAM 

Without  delving  more  into  the  theory  itself,  which  can  be 
found  in  Ref.  16,  let  us  apply  it  to  the  same  line  source  as  in 
the  previous  section.  The  transmission-image  source  can  be 
written  in  the  form 

J.O,  p)  -  6,/fMp)  +  F{(Bh’,  p)|«(y  +  jb  sin  8)S(z  -  hr), 

(28) 


where  the  function  F\'  is  defined  as 

F,'<Bh',  p)  -  H{y~h,}  [(H  +  fc')V,(p)  +  (H  -  h')V3(p)], 

(29) 


and 


H  *  H(h',  p)  «  [h'z  +  (p/B)2\in, 

h'  ■  h+  jb  cot  8.  (30) 

The  convergence  condition  Im[p]  ”  0  defines  a  path  z  » 
H(p),  which  consists  of  a  part  of  a  hyperbolic  curve  in  the 
complex  z  plane.16  To  obtain  an  asymptotic  result  from  the 
exact-image  field  solution,  let  us  once  again  study  the  far 
field  of  the  beam.  After  making  use  of  expression  (12),  we 
must  now  apply  the  following  approximation  for  large  IH, 
instead  of  expressions  (13),  for  x  *  0,  z  <  0: 

D  *  D0  +  qH{p),  D0  -  I0v  +  Jb  »in  9)2  +  z2]1'2, 

q  ■  lzl/D0  *  cos  S2.  (31) 

Thus  the  field  in  medium  2  can  be  written  in  the  approxi¬ 
mate  form 

£0)  *  -  ju^UJG^Do)  jl  +  F,'(Bh\  p) 

x  exp(-;kjfl//(p)]dpl  .  (32) 
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Fig.  3.  Gaussian  beam  transmitted  through  the  interface  seen  as 
arising  from  an  approximate  point  source  at  point  uh2  +  jS2. 


Here,  the  Green  function  is  defined  as  for  Eqs.  (S),  for  the 
medium  2  with  k2  replacing  k{.  Applying  an  integral  identi¬ 
ty  (Ref.  16,  Eq.  (17)|,  we  can  write,  for  the  expression  in 
braces  in  expression  (32), 

1  +  [  F \'(Bh',  p)  exp(-yy/f(p)]dp  -  — ■  2 7 

Jo  y  +  (7*  -  r  * 

npH(y*  -  B2h')l/iJ.  (33) 

Use  of  this  and  (see  Fig.  3  for  the  different  parameters) 

D0  *  r  +  jb  sin  0  sin  02  »  r2  +  d'  sin  02,  (34) 

with 

d'  »  d  +  jb  sin  0,  (35) 

makes  it  possible  to  write  expression  (32)  in  the  following 
form: 


1  ,  ,/  2;  \w*  2k2  cos  02 


X  exp «xp(-/k,  sin  0,d')  exp (~jki  cos  0,h').  (36) 


Here,  we  have  made  the  far-fiehl  assumption  r2  »  h,  d,  b; 
otherwise  the  expressions  would  have  been  much  more  com¬ 
plicated.  Expression  (36)  can  be  interpreted  in  terms  of 
geometrical-optics  rays  of  the  Gaussian  beam,  with  a  phase 
shift  in  each  medium  and  a  transmission  coefficient  at  the 
interface.  From  this,  an  equivalent-source-point  location 
can  be  defined  from  the  divergence  of  two  adjacent  rays  close 
to  the  beam  axis  9\  ■  0.  Taking  the  combination  of  the  last 
two  exponentials  in  expression  (36),  we  may  require  that  the 
corresponding  expression  for  the  image  case  have  the  same 
change  when  the  angle  0|  differs  from  the  axial  angle  9. 
Writing  for  the  exponential  expressions 


h  h 

+1  =  -jkx(d'  sin  0,  +  h'  cos  0^)  *  -j  — —  +  *  cos(0,  -  0), 

aas D  1  1 


(37) 

*2  *  ->k2(d2'  sin  02  +  h,'  cos  0,) 

k.j/i2 

"  -J  C08  f  +  *2*2  COS(02  0O), 

(38) 

k,  sin  0  ■  k2  sin  90, 

(39) 

we  may  require  that  phase  and  amplitude  changes  along  the 
plane  z  ■  0  be  the  same  for  the  original  case  and  the  image 
case.  This  implies  that  the  differentials  of  * ,  and  +2  are  the 
same: 

d*!  *  d*2. 

(40) 

When  supplemented  by  the  condition  between  the  differen¬ 
tials  d£t  and  d02  arising  from  Snell's  law,  A,  sin  9\  -  k2  sin  02, 
Eq.  (40)  gives  us  the  following  expression  for  the  unknown 

hi¬ 


lt ,  cos3  0, 

h2«h- - — ^ 

k,  cos3  0, 


at  0,  ■  0. 


(41) 


Because  of  L\  »  hi cos  0i,  L2  “  hj/co»  02,  the  corresponding 
result  for  L2  is 


2cos202 
k,  cos2  02 


(42) 


which  was  also  given  in  Ref.  5.  To  obtain  an  expression  for 
62,  defining  the  apparent-source  location  in  complex  space, 
we  must  study  the  real  parts  of  Eqs.  (37)  and  (38).  The 
direction  of  the  vector  is  known  from  the  main  beam  direc¬ 
tion  in  the  half-space  2.  The  magnitude  62  is  not  obtained 
from  Eq.  (40),  since  the  first  differentials  vanish  identically 
at  0t  «  0.  Taking  the  second  differentials  of  *1  and  ♦2  and 
equating  the  real  parts  gives  us  the  result 

k,  cos2  02 

b,  -  6  — - —=■  for  0,  «  0,  (43) 

k  |  cos2  0i 


a  condition  that  was  also  given  in  Ref.  5. 

Thus  the  location  of  the  approximate-image  source  can  be 
obtained  through  two  different  approaches:  the  present 
exact-image  theory  and  the  saddle-point  asymptotic  analy¬ 
sis  of  the  Fourier-integral  representation  of  the  transmitted 
field  in  Ref.  5.  The  present  approach  has  the  advantage  of 
working  with  sources,  which  gives  a  more-physical  insight 
for  the  problem,  in  addition  to  being  exact. 


CONCLUSION 

In  this  paper,  the  exact-image  method,  which  was  previously 
introduced  for  problems  involving  sources  in  real  space  with 
two  homogeneous  media  and  a  planar  interface,  is  extended 
to  problems  involving  sources  in  complex  space.  This  is 
possible  because  the  analytical  functions  of  real-space 
source  poeition  can  be  extended  to  functions  of  complex- 
space  source  position,  an  idea  that  has  been  applied  before 
with  success  in  diffraction  problems.  The  complex-source- 
point  theory  makes  it  possible  to  analyze  Gaussian -beam 
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problems  with  a  simple  source  instead  of  a  complicated 
source  in  real  space.  The  method  has  been  tested  in  this 
paper  by  analyzing  asymptotic  (far-field)  expressions  for  a 
reflected  Gaussian  beam  and  a  transmitted  one,  resulting  in 
well-known  Goos-Hanchen  and  angular  shifts  for  the  re¬ 
flected  beam  and  in  an  apparent  position  for  the  transmitted 
beam.  The  method  described  here  is  recommended  for 
more-exact  calculation  of  the  fields  in  these  beams,  for  ex¬ 
ample,  when  the  asymptotic  approximations  will  not  work. 
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Theory  of  Time-Domain  Quasi-TEM  Modes 
in  Inhomogeneous  Multiconductor  Lines 

ISMO  V.  LINDELL.  senior  member,  ieee,  and  Q1ZHENG  GU 


—Theory  for  quasi-TEM  modes  propagating  in  a  transversely 
■imuumi—  (muitidMectnc),  longitudinafiy  uniform  transmission  line. 


derived  for  time- harmonic  eaves,  is  derived  for  transient  sig- 
J***  u  men  that,  while  the  starting  point  for  the  theory  is  completely 
the  result  is  similar  to  the  rime- harmonic  theory,  and  previously 
fnftmn  for  propagating  modes  also  apply  in  the  transient  case, 
of  applicability  is  discussed  with  a  simple  example. 

I.  Introduction 

THE  INHOMOGENEOUS  multiconductor  transmis¬ 
sion  line  is  interesting  because  of  the  wide  field  of 
plications  0f  the  microstrip  structure.  In  the  case  of 
nie-harmonic  (sinusoidal)  fields,  the  exact  TEM  modes 
for  this  kind  of  geometry,  do  not  exist.  Before  the  intro¬ 
duction  of  a  theory  for  the  quasi-TF.M  modes,  there 
were  often  erroneous  assumptions  about  the  nature  of 
fields  propagating  in  such  structures.  The  exact  theory 
^  two-conductor  lines  was  first  given  by  dos  Santos  and 
Figaiuer  in  1975  [1],  generalized  to  multiconductor  lines 
by  Mannersalo  in  1977  (not  published)  and.  in  a  more 
complete  form,  by  Lindell  in  1981  (2).  There  also  exist 
other  expositions  on  the  subject  (3),  [4). 

As  pointed  out  in  [2],  the  transverse  quasi-TEM  mode 
field  is  a  kind  of  glued-together  pair  of  static  electric  and 
static  magnetic  fields,  whose  boundary  conditions  are  con¬ 
nected  through  a  pair  of  consistency  equations,  which 
also  determine  the  propagation  properties  of  the  quasi- 
TEM  wave.  The  longitudinal  components  can  be  calcu¬ 
lated  from  the  transverse  components  and  they  are  low  for 
small  frequencies.  This  theory  was  obtained  by  expanding 
all  unknown  quantities  in  series  of  ascending  powers  of  the 
frequency.  How  this  can  be  generalized  to  fields  of  arbi¬ 
trary  time  dependence  is  not  evident,  because  there  does 
not  exist  a  small  parameter  like  w.  However,  obviously,  a 
umilar  theory  should  apply  for  transient  signals  whose 
vpectrum  consists  of  sufficiently  low  frequencies.  A  need 
for  a  firm  foundation  for  such  a  theory  exists  because  of 
the  new  generation  of  computers  applying  microstriplike 
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geometry  and  rapid  pulsed  signals,  although  the  important 
wavelengths  of  such  signals  are  still  large  with  respect  to 
transverse  measures  of  the  line.  Calculations  of  coupling 
between  such  signals  applying  the  ideal  TEM  mode  ap¬ 
proximation  have  recently  been  presented  [5),  [6]. 

The  theory  given  here  is  not  based  on  an  asymptotic 
(perturbational)  series  approach,  but  on  an  iterational 
consideration,  which  starts  on  a  quasi-TEM  assumption, 
proceeds  to  derive  properties  for  such  a  field,  and  finally 
tries  to  find  out  conditions  of  validity  for  the  original 
assumption. 

II.  The  Iteration  Method 

The  following  theory  of  mode  propagation  in  a  trans¬ 
versely  inhomogeneous  multiconductor  waveguide.  Fig.  1, 
is  a  companion  to  that  given  in  [2]  and  uses  much  of  the 
same  notation,  the  main  difference  being  that  the  sub¬ 
scripts  z  are  deleted.  The  electromagnetic  field  is  studied 
in  terms  of  its  transverse  and  longitudinal  components, 
respectively,  transverse  and  parallel  to  the  guide  axis  vec¬ 
tor  u: 

E(r,t)  **e{r,t)+ue(r,t)  (1) 

H(r,t)  -h(r,t)+uh{r,t).  (2) 

When  the  Maxwell  equations  are  written  for  the  decom¬ 
posed  fields  and  the  resulting  equations  are  decomposed 
into  axial  and  transverse  components,  we  have  the  set  of 
equations 


Vj.  Xr»  -  und,h 

(3) 

V  ±  X  h  =  urd,e 

(4) 

V±e  =  -  d.e  +  d,u  x(pA) 

(5) 

V ±h  -  -  d.h  -  d,u  X(d) 

(6) 

Vj_  (ce)  -  -  id.e 

(7) 

Vx  (nh)  -  -n9th. 

(8) 

Here,  vA  denotes  the  transverse  component  of  v;  p  and  e 
are  functions  of  the  transverse  vector  p,  and  the  partial 
derivatives  dF/d f  are  denoted  in  short  by  d(F. 

Equations  (3)-(8)  are  exact.  As  a  starting  point  to  the 
iteration  method  leading  to  the  quasi-TEM  theory,  we 
assume  that  the  right-hand  sides  of  (3),  (4),  (7),  and  (8)  are 
small  and  approximate  them  by  zeros.  This  assumption 
will  be  studied  more  closely  later  on.  After  solving  for  the 
transverse  fields  on  the  left-hand  side,  the  longitudinal 
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Fig.  1.  Cross  section  of  the  transversely  inhomogeneous  multiconductor 
line  with  N  conductors  S,  and  a  conducting  sheath  $>.  The  unit  vector 
5  and  t  coordinate  are  perpendicular  to  the  plane. 


fields  are  solved  from  (5)  and  (6).  The  iteration  could  be 
continued  by  substituting  these  on  the  right-hand  sides  of 
(3),  (4),  (7),  and  (8),  and  solving  for  the  next  approxima¬ 
tions  for  the  transverse  fields.  The  process  will  evidently 
converge  if  the  longitudinal  fields  are  small  enough.  Here, 
we  are  only  concerned  with  the  first  round  of  iteration, 
which  gives  us  static  approximation  of  the  transverse  fields 
and  quasi-static  approximation  of  longitudinal  fields. 

III.  Static  Approximation 

Setting  e-0  and  A-0  in  the  four  equations  (3),  (4), 
(7),  and  (8)  gives  us  equations  for  the  first  approximation 
of  the  transverse  fields  e0,  A0; 

Vxxio(r,t)=0  (9) 

Vx(te0)=0  (10) 

Vxxh0(r,t)=  0  (11) 

Vo.-(mSo)~0.  (12) 

Equations  (9)  and  (10)  constitute  an  electrostatic  problem 
in  two  dimensions,  and  as  the  boundary  conditions  on  the 
conductors  we  have 

nxe0«0.  (13) 

Correspondingly,  (11)  and  (12)  define  a  magnetostatic 
problem  with  the  boundary  conditions 

n-Ao-0.  (14) 

Here,  n  denotes  the  unit  vector  normal  to  the  conductor 
boundaries.  The  static  problems  can  be  treated  with  poten¬ 
tial  quantities.  In  fact,  defining  the  electric  field  in  terms 
of  a  scalar  potential  <p: 

e0(r,t)  -  -  V^(r,f)  (15) 
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•  °CT0*»  1^ 

means  the  voltage  between  the  ith  conductor  un 
sheath.  The  present  theory  is  not  limited  to  closed*""  *** 
tries.  If  the  correct  behavior  of  potential  functions8”"!!! 
infinity  is  taken  into  account,  the  analysis  can  be  ^ 
in  very  much  the  same  fashion.  However,  instead  oTtfeT 
ning  voltages  with  respect  to  the  sheath,  another  rtfttaL 
potential  (“ground”)  must  be  assumed. 

The  solution  of  (16)  and  (17)  obviously  depends  on  the 
set  of  boundary  values  U,(z,i),  which  can  be  represented 
as  an  N-vector  U(z,  /).  To  be  able  to  write  the  solution  for 
any  boundary  values,  we  must  solve  N  different  normal- 
ized  problems  for  functions  <p,(p),  making  up  an  \-veaor 
function  4>(p),  satisfying  the  boundary  conditions  +  -i 
on  Sj  and  -  0  on  all  the  other  conductors  (i  *  j).  The 
solution  corresponding  to  the  boundary  value  vector  U(:.i) 
can  then  be  written  as 
v 

E  VAz,t)$,(p)  ~U(z, f )•  ♦(p).  (1*1 

i-i 

Here,  we  have  adopted  the  notation  °  for  the  inner 
product  of  two  Af-ve ctors  to  distinguish  it  from  the  “dot* 
product  of  two  vectors  in  the  physical  three-space. 

In  the  same  manner,  the  magnetostatic  problem  can  be 
formulated  in  terms  of  an  axial  vector  potential  A  »  uA: 

p(p)h0(r,t)  =Vxd(r,l)xu.  (W) 

Equation  (12)  is  automatically  satisfied  through  (19).  Out 
of  (11)  and  (14)  comes 

<20) 


(9)  will  be  automatically  satisfied,  whence  the  equations 
for  <(>,  from  (10)  and  (13),  are 

VA(*y»-0  (16) 

<p(r,t)  for peS,.  (17) 

Here.  S,  denotes  the  surface  of  the  ith  conductor,  i  - 
0, 1,  •  ■  • ,  N.  For  simplicity,  we  assume  that  there  is  a  sheath 
conductor  denoted  by  i »  0  and  define  U0  »  0  so  that  U, 


A(r,  t)  * ♦,(*, /)  forpeS,. 

Here,  denotes  the  magnetic  flux/unit  length  between 
the  conductor  i  and  the  sheath.  In  terms  of  a  genet* 
solution  A'-vector  A(p)  and  a  set  of  boundary  value* 
♦(z,  /)  the  solution  can  be  written  as 

A(r.t)  -ir(z,t)»  A(p).  (22) 

In  the  electrostatic  problem  (16),  (17),  the  bound*7 
value  voltage  Ai-vector  uniquely  defines  the  quantities 
i.e.,  the  charge/unit  length  on  each  conductor  through 

Q,(z.  t )  =<^  te0  ndl  *  ~^> «—  dl  ^ 

where  C,  denotes  the  circumference  of  the  i  th  con<*“c^ 
The  linear  relation  between  the  voltage  and  charge  ‘  . 
tors  can  be  written  in  terms  of  a  static  capacitance/ 
length  matrix  C 

( 2 4) 

Q(z.t)  •C»(J{z.t). 

For  the  magnetostatic  problem,  there  is  also  a 
relation  between  the  magnetic  fluxes  and  the  currents  , 
the  conductors,  defined  by 

t-  -  t  1  iA  M  (25) 

j,-r^ 

tf(r,f)«L<>/(z.r) 

where  L  is  the  inductance/unit  length  matrix. 


IV.  The  Quasi-TEM  Field 


Until  now,  the  two  static  problems  have  no  connection. 
■^e  boundary  value  .V-vectors  are.  however,  coupled 
through  (5)  and  (6).  From  (5).  we  can  solve  the  approxi- 
^le  longitudinal  component  e,: 

el(r,t)~3.^(r,t)+d,A(r,t).  (27) 


corresponding  approximation  for  the  longitudinal 
magn«tic  field  A,  cannot  be  obtained  from  (6)  explicitly: 

(r.t)  =  -  3.h0(r.t)-  u  xtd,e0(?,t) 


•  u  X 


1 

—  Vx  d.A(r,  / )  +  «V  x  d,<b(  r,  t ) 


•  (28) 


Because  the  right-hand  side  of  (28)  is  curl  free,  as  can 
easily  be  checked.  A,  can  be  obtained  through  integration 
from  a  reference  point  p0  to  the  general  point  p.  Integrat¬ 
ing  round  the  conductor  i,  the  integral  should  vanish  since 
h  iS  a  unique  physical  quantity.  To  obtain  unique  values 
for  A,  from  (28)  by  integration,  one  more  condition  for  hx 
is  needed,  because  the  reference  value  At(  p0)  is  not  known. 
The  missing  condition  is  obtained  by  integrating  (3)  over 
[he  waveguide  cross  section,  which  gives  us  zero  if  Stokes’ 
theorem  is  invoked,  because  the  line  integral  of  the  electric 
field  around  the  perimeter  vanishes.  Thus,  d,jhdS  -  0,  or 
the  integral  is  constant  in  time.  If  the  fields  are  switched 
on  at  some  finite  time,  the  constant  is  zero  and  we  have  for 
the  additional  condition  for  A, 

fhl(r.t)dS~0  (29) 

Js 

where  S  is  the  total  cross  section  surface  of  the  waveguide. 

Because  of  the  boundary  condition  e, »  0  on  the  con¬ 
ductors,  (27)  defines  a  relation  between  the  boundary 
values  of  the  potentials  <p  and  A : 


d,f/(z./)+d,*(z.r)=0.  (30) 


This  is,  in  fact,  one  form  of  writing  the  Faraday  law,  in 
which  a  moving  magnetic  flux  produces  an  induced  volt¬ 
age.  The  other  equation,  (28),  gives  rise  to  another  relation 
between  integrated  quantities  of  the  potentials  <f>  and  A,  in 
(he  form  of  a  continuity  equation: 

a./(z.r)  +  d, <?(*,/)-  0.  (31) 

Substituting  from  (24)  and  (26),  we  have 

(32) 

*./(*,/)--£•  3, t/(z.r).  (33) 

T'1es*  are  recognized  as  transmission-line  equations  for  a 
"•nlticonductor  line.  Eliminating  /,  we  have  the  wave 
aquation 


{dH~  k°Cd^)oU{2.t)=  0.  (34) 

y 

d?n°‘«  t*,e  iV  *  'v  un‘l  matrix-  Limiting  ourselves 
so  ution  traveling  in  the  positive  :  direction,  the 
°r ,n  ( 34)  can  be  halved  to  produce  the  equation 

[dJ+(L'Q)'/2d,)'U{:.t)-0.  (35) 


Because  C  and  L  are  positive  definite  and  symmetric 
matrices,  there  exists  a  square-root  matrix  ( L  •  C)l/J  which 
is  positive  definite.  Hence,  it  possesses  a  complete  set  of 
eigenvectors  and  positive  eigenvalues,  and  the  solution  of 
(35)  can  be  written  in  terms  of  the  solutions  of 

(L-C)1/2»F'  =  -F'.  (36) 

Corresponding  to  each  eigenvector  V>  there  exists  a 
mode  of  the  original  problem,  whose  voltage  vector  can  be 
written  as 

UJ(z,t)  -V‘f\z  -  Ujt)  (37) 

with  an  arbitrary  function  fr  Thus,  the  eigenvalues  of  the 
algebraic  equation  (36)  define  the  velocities  of  the  S 
quasi-TEM  modes  on  the  inhomogeneous  waveguide.  The 
current  distribution  vector  corresponding  to  the  y  th  mode, 
1‘.  satisfies  similar  equations  but  with  L  and  C  inter¬ 
changed.  As  in  [2],  this  leads  to  an  eigenvector  different 
from  U1,  in  general,  which  means  that  there  does  not  exist 
a  simple  scalar  impedance,  but  the  impedance  is  a  matrix 
quantity,  denoted  here  by  Z. 

V.  Traveling  Wave  Solutions 

The  practical  multiconductor  line  problem,  with  possible 
small  deviation  from  axial  uniformity,  can  be  treated  with 
traveling  wave  quantities  with  less  effort  than  with  cur¬ 
rent/voltage  quantities.  In  fact,  we  can  define  the  voltage 
wave  vectors  (7)  by 

Ut-U±Z°!  (38) 

where  the  impedance  matrix  is  defined  in  any  of  the 
following  equivalent  forms  [2]: 

*(40£)  l/3‘,4a“4°(£0^)  l/2-  (39) 

It  must  be  noted  that  since  in  the  general  case  the  matrices 
C  and  L  do  not  commute,  the  square  root  of  their  product 
is  dependent  on  the  order  of  the  product.  Solving  the 
voltage  and  current  vectors  from  (38)  in  terms  of  voltage 
waves  and  substituting  in  (32)  give  us  the  two  equations 

d]Ut  ( z,i)±(L  •  C)*/2  •  d,Ut  (z,  t )  -  0.  (40) 

It  is  seen  that  equations  (40)  are  uncoupled  for  both 
traveling  voltage  vectors.  Thus,  they  possess  solutions  of 
the  type 

f/'(r,r)«F'/,(zTV).  (41) 

For  a  line  with  slight  nonuniformity  along  its  axis, 
equations  (40)  can  be  generalized  to  possess  a  perturba- 
tional  coupling  term  on  the  right-hand  side  (5),  (6). 

It  is  also  easy  to  show  that  the  integral  of  the  Poynting 
vector  ^e0  x  A  J  over  the  cross  section  equals  the  sum  of 
:  U  •  /*.  or  the  propagating  power  in  the  quasi-TEM  mode 
can  be  expressed  in  terms  of  fields  or  boundary  values. 
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VI.  On  the  Validity  of  the  Quasi-TEM 
Concept 

After  having  calculated  the  approximate  longitudinal 
field  components  in  terms  of  the  approximate  transverse 
components,  we  can  substitute  them  on  the  right-hand 
sides  of  (3),  (4),  (7),  and  (8)  and  calculate^  a  better  ap¬ 
proximation  for  the  transverse  ftelds,  e2,  A2.  For  these 
fields  not  to  differ  from  the  original  approximations  e0,  h0 
significantly,  the  right-hand  sides  should  be  small.  What 
does  this  mean?  Obviously,  in  order  to  be  able  to  ap¬ 
proximate  the  transverse  curl  operations  in  (3)  and  (4)  by 
zero,  the  right-hand  sides  should  be  small  with  respect  to 
other  derivatives  than  curl  of  the  vector  functions  on  the 
left-hand  sides.  In  this  case,  the  iteration  would  converge 
and  the  terms  already  calculated  would  present  a  reason¬ 
able  approximation  to  the  fields. 

Let  us  study  the  question  in  terms  of  simple  considera¬ 
tion  of  order.  If.  for  example,  the  following  inequality  is 
valid: 

!IVa«0||»|m«,*iI  (42) 

then  obviously  from  (3),  V  x  e0  is  small  when  compared 
to  other  transverse  derivatives  of  e0  and  it  can  be  ap¬ 
proximated  by  zero.  Thus,  when  solving  the  next  itera- 
tional  step,  the  transverse  field  e2  in  terms  of  A,  from  (3). 
the  solution  e2  does  not  differ  from  e0  very  much.  Now 
the  order  of  Vj.  is  1/D,  where  D  is  the  transverse 
dimension  of  the  guide.  If  the  dimension  of  the  axial 
variation  of  the  fields  is  L.  we  have  /’( z  -  vpt)  of  the  order 
f(z  -  vpt)/L.  where  vp  is  the  phase  velocity  of  the  mode 
and  from  (28)  and  (42)  we  may  write 

l*ol  i-*o  + (43) 

This  inequality  is  obviously  valid  if  D/L  is  small  enough, 
i.e..  if  the  transverse  dimension  of  the  guide  is  small 
enough  with  respect  to  field  variations  in  the  z  direction. 
In  other  words,  the  variation  of  the  signal  should  be  so 
slow  that  the  propagating  field  does  not  change  much  at 
the  distance  of  the  transverse  dimension  in  z  direction.  But 
this  is  not  all.  because  the  vector  term  on  the  right-hand 
side  may  also  be  small.  This  happens  if  the  inhomogeneity 
of  the  guide  is  small.  In  fact,  if  the  guide  becomes  homoge¬ 
neous.  the  solution  of  (20)  becomes  a  multiple  of  the 
solution  of  (16)  and  the  relation  obtained  from  (30)  is 
$  -  vpA .  or  from  (15)  and  (19)  h0  -  u  x  e0/tj,  with  ij  » 
pvp  « l/«up,  which  makes  the  right-hand  side  of  (43)  vanish. 
If  the  inhomogeneity  is  small,  the  term  need  not  vanish, 
but  it  is  small  and  (43)  is  valid  for  larger  values  of  the  ratio 
D/L. 

There  does  not  seem  to  be  an  easy  way  to  express  the 
condition  of  validity  in  more  exact  terms.  A  simpler  condi¬ 
tion  is  obtained  if  the  ratios  |eii/|e0|,  |A,|/A0|  are  consid¬ 
ered.  This  not  only  shows  us  that  the  basic  approximation 
of  the  field  is  quasi  TEM.  but  also  gives  an  idea  on  the 


octomrihj 


Fig.  2.  The  inhomogeneous  coexist  cable.  Basic  wave  is  a  TE  vase, 
which  is  quasi-TEM  for  sufficiently  slowly  varying  transieni  sigaab 
and/or  small  inhomogeneity  (small  t,  - 1). 


convergence.  In  fact,  writing 

k,|  _  f'(z-vpt)  |4» -  o„A\  _  f  I(4>-'*)°F1 
|e0l  "  f(z~vpt)  IV  >|  “  /  |  V±fV\ 

we  may  require  this  to  satisfy  the  condition  •«  1.  The  last 
factor  is  zero  for  homogeneous  line  and  small  for  small 
inhomogeneity.  Because  the  maximum  values  of  4  and  A 
are  their  boundary  value  1,  it  is  obviously  bounded  by  the 
value  2D,  which  shows  us  that  the  quasi-TEM  mode  is 
always  possible  if  the  signal  is  varying  slowly  enough. 


VII.  Example 

Let  us  elucidate  the  problem  of  convergence  through  a 
simple  example  of  inhomogeneous  coaxial  cable  with  p  m  f*o 
and  i,  - 1  for  8  -  0  ■  ■  ■  ir  and  *,  *  1  for  8  «*  0  •  •  •  _  *•  | 

2.  Because  of  the  special  symmetry,  the  potential  fiehk  ♦ 
and  A  are  multiples  of  the  same  function  g(p): 


lb\ 

In  —  I 

\  P  I 


The  inductance  is  the  same  as  for  the  homogeneous  < 
line: 


L-Ln-^ln  - 


and  the  capacitance  equals 


l+<,\  ir(l  +  «,)<0 

‘r,C°'  '»(-)  ' 
a  I 


The  phase  velocity  of  the  quasi-TEM  wave  is 


r 
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The  potentials  are 

<t>(r,t)  =  Vg(p  )f(z-  vpt)  (49) 

A(f,t)~'*gil>)f(z-  V)  "  —g(p)f{z-vpt).  (50) 

UP 

from  (27)  we  have 

V 

l{,Vg(p)f{z- v?t)~  v—g{p)f(z  -  <y)  =0  (51) 

the  longitudinal  electric  field  vanishes  in  this  approxi- 
u^tion.  Thus,  the  quasi-TEM  is  in  fact  a  TE  field  in  this 
proximation,  which  is  due  to  the  special  symmetry  of  the 
juncture.  From  (28)  we  have 

k(r.<)-”/'(*- V)/  f ~ —  «r(p)  —  ) 

n  J>0\vp  c  i 

(fn-V,g)dC  +  hi(r0).  (52) 

unit  vector  m  is  normal  to  the  integration  path. 
Qcrl...tf  m-Vj.gip)  =  0  on  each  radial  lin  ,  it  is  clear  that 
^  jS  a  function  of  9  only.  Thus,  the  integration  path  can 
^  taken  along  any  circle  of  radius  p.  If  the  constant  h,( r0) 
n  determined  through  the  condition  (29),  the  result  is 

for  -  ir  «  9  <  it.  (53) 

This  expression  is  seen  to  vanish  for «,  »  1,  in  which  case 
the  guide  is  homogeneous.  The  null  field  is  obtained  for 
|»i*/2  and  maxima  at  9  =  0. ir.  To  check  the  quasi- 
TEM  character  of  the  wave,  we  compare  the  magnitudes  of 
*,  and  h0: 

|A,|  /'  e,  —  1  \  f'(z  -  v  t)  I  ir  \ 

f(z~  V)  |  i]0]~  2)4  (54) 

The  maximum  value  of  the  last  factor  is  irb/ 2.  Thus,  the 
right-hand  side  of  (54)  gives  us  the  relative  rate  of  change 
of  the  signal  in  the  axial  direction  over  the  distance 
i«,-l)ir6/2(«,  +  1).  Obviously,  the  convergence  of  the 
iteration  is  good  and  the  quasi-TEM  concept  valid  if  this 
rate  is  small. 

VIII.  Conclusions 

The  quasi-TEM  mode  theory  of  inhomogeneous  multi- 
“nductor  waveguides,  previously  presented  for  time- 
^'"’onic  Helds,  was  generalized  for  fields  with  arbitrary 
dependence.  The  theory  is  based  upon  an  iterative 
|  and  the  condition  for  its  convergence  was  out- 

and  elucidated  with  a  simple  example.  It  was  seen 
d*  resu,t‘nS  ‘heory,  basically  similar  to  that  given  for 
harmonic  case,  can  be  applied  for  transient  sig- 

iht  ^  var*a,'on  of  the  signal  is  slow  enough  or 

"Homogeneity  is  not  too  large. 
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Abstract— Tht  aaalysh  of  rooMnce,  input  impedance,  and  rudixikm  of 
the  cWptk  disk  nkrostrip  structure  b  rigorously  formulated,  using  the 
scalar  and  vector  Malhico  transforms.  With  the  help  of  these  transforms, 
the  resonance  frequencies  of  the  structure  can  be  derived  exactly  using 
GaierUa’s  method  and  approximately  using  a  perturbationai  approach. 
Expressions  for  the  input  impedance  and  the  radiation  pattern  are  also 
obtained. 


I.  Introduction 

IT  IS  KNOWN  that  circular  or  rectangular  disk  microstrip 
antennas  usually  radiate  waves  which  are  linearly  polarized 
[l]-[3].  However,  in  such  structures  circular  polarization  can 
be  obtained  using  a  complicated  feed  setup  comprising  more 
than  one  teed. 

From  experimental  measurements  (I),  [2|  and  recent 
analytical  work  [3],  it  is  shown  that  by  using  a  slightly 
elliptical  shaped  disk,  circular  polarization  can  be  achieved 
while  retaining  a  single,  simple  feed.  Also,  for  applications  in 
harmonic  multipliers  and  parametric  amplifiers,  the  circular 
microstrip  disk  resonator  is  unsuitable  due  to  a  harmonious 
relationship  of  mode  frequencies.  However,  suitable  modes 
for  both  applications  car  be  found  in  the  resonances  of  an 
elliptic  structure  where  a  further  degree  of  freedom,  the 
eccentricity,  enables  the  desired  frequency  relationships  to  be 
achieved  |4)-£6J.  Thus,  the  eccentricity  as  a  design  parameter 
provides  an  additional  flexibility  and  enhances  the  usefulness 
of  the  elliptic  disk  structure.  An  added  advantage  of  the  use  of 
the  elliptic  structure  over  the  circular  one  is  that  the  former  can 
be  designed  to  have  a  wider  frequency  band  of  operation  than 
the  latter.  A  broad-band  structure  can  be  realized  out  of  the 
circular  structure  only  through  the  use  of  more  than  one 
element  [7H9J. 

In  past  studies  of  the  elliptic  disk  microstrip  structure,  more 
work  was  devoted  to  the  analysis  of  the  structure  as  a  resonator 
(4]-[6]  than  as  an  antenna  [i]-{3].  In  [3],  the  analysis 
of  the  elliptic  resonator  is  carried  out  by  employing 
the  magnetic  wall  cavity  model,  whereas  in  (6),  the  analysis  is 
carried  out  using  a  quasistatic  approximation  where  an 
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equivalent  capacitance  of  the  elliptic  structure  is  derived.  In 
both  analyses,  the  resonant  frequencies  obtained  are  pure  real 
and  therefore  do  not  account  for  the  radiation  losses,  thus 
rendering  the  obtained  results  very  limited.  Also,  using  the 
equivalent  capacitance  to  calculate  the  resonant  frequencies  is 
proven  to  be  incorrect  [  10]— [12]. 

In  the  analysis  of  the  elliptic  structure  as  a  radiator  and  in 
calculating  the  radiation  pattern  [3],  the  fringing  effects  of  the 
electric  field  are  taken  into  account  by  superficially  imposing 
an  impedance  boundary  condition  at  the  opening  of  the  cavity. 
This  surface  impedance  is  approximated  by  the  one  obtained 
from  solving  the  circular  disk  antenna,  thus  limiting  the 
obtained  results  for  small  eccentricity. 

In  this  paper,  a  rigorous  analysis  of  the  elliptic  microstrip 
structure  is  carried  out  using  the  scalar  and  vector  Mathieu 
transforms  [13].  These  transforms  allow  an  exact  formulation 
of  the  elliptic  disk  structure  which  has  long  been  thought 
almost  impossible  to  formulate  rigorously  [6].  h  is  shown  that 
the  current  distribution  on  the  disk  is  rigorously  derived  using 
these  transforms  and  that  it  is  governed  by  a  set  of  vector 
integral  equations.  These  equations  are  then  solved  using  the 
Galerkin's  method,  in  which  the  current  is  expanded  in  terms 
of  the  orthogonal  set  of  the  transverse  magnetic  (TM)  and 
transverse  electric  (TE)  modes  of  the  magnetic  wall  cavity, 
since  they  form  a  complete  set  of  basis  functions,  to  obtain 
linear  algebraic  equations  in  the  expansion  coefficients.  The 
eigenvalue  equation  is  then  obtained  by  setting  the  determinant 
of  the  coefficients  of  these  algebraic  equations  to  zero.  In  the 
limit  of  a  thin  substrate,  a  perturbative  approach  is  used  to  get 
a  simple  expression  for  the  resonant  frequencies.  Finally,  the 
input  impedance  and  the  radiation  pattern  are  derived  both 
exactly  and  in  the  limit  of  a  thin  substrate. 

II.  Field  Expressions  Due  to  the  Current  Distribution  on 
the  Disk 

Fig.  1  shows  the  geometrical  configuration  of  the  problem. 
A  perfectly  conducting  elliptical  disk  is  placed  on  top  of  a 
dielectric  substrate  backed  by  a  perfectly  conducting  ground 
plane.  The  disk  has  a  semimajor  axis  a  and  a  semiminor  axis 
b  and  of  a  focal  length  2c0.  The  dielectric  substrate  has  a 
permittivity  of  «,  and  thickness  d  An  elliptical  coordinate 
system  (w,  it,  z)  is  used  to  express  points  in  space.  These  are 
related  to  the  Cartesian  coordinates  as  follows: 

jr=Co  cosh  (m)  cos  (t>),  y  =  c0  sinh  («)  sin  (v). 

The  disk  is  fed  by  a  probe  P  situated  at  (u0,  u0).  Using  the 
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stratified  medium  formulation  [  10]-[  12],  (14J,  the  z-compo- 
nent  of  the  electric  and  magnetic  fields  in  the  free  space  region 
due  to  the  current  distribution  on  the  disk  can  be  shown  to  be 
given  by 


£*=E  j0  ****** 

fi.a 


-  R™e‘ktt2H**,]}[iot„(c0kl„  u,  t>)  (1) 
=  S  j0  dk'k,  ha„(k,){et,ki: 

m,a 

+  Rne'ktaH+t'WctAeokf,  u,  u)  (2) 

where  the  positive  sign  is  for  z  >  0  and  the  negative  sign  is  for 
z  <  0.  \l>a„(c0k„  u,  u)  is  an  elliptic  harmonic  given  by  the 
product  of  a  radial  and  an  angular  Mathieu  function  and  is 
defined  in  [13,  appendix  A],  H  is  the  height  of  the  disk  above 
the  substrate  which  will  be  later  set  equal  to  zero  and  is 
introduced  to  avoid  the  ambiguity  in  applying  the  boundary 
conditions  at  the  plane  of  the  disk.  The  transverse  components 
of  the  electric  and  magnetic  fields  can  be  readily  obtained  from 
Ez  and  Hz  using  the  following  formulas  [14]: 

£,(*,,  u,  v,  Z)  =  ^i  j^r  ( k„u ,  u,  z)^ 

+  i<jiHoVTx(ifIz(kt,u,  u,  z))j 

fM*,,  u,  v,  z)  =  —  £vr  ( k,,u ,  v.  z)J 

-iu«Vrx(i£s(ke,u,  v,  z„] 


and  toa^Cok,,  u,  v)  is  given  by  [13] 
v) 


d  d 

1  iu  *a^c°kl’'  ^  +°*<cok"  u<  t'> 

h  d  a 

—  \ka„(Cok„,  u,  v)  \l>a„(c0k,,,  u,  v) 

av  du 


The  transverse  electric  field  due  to  the  current  distribution 
on  the  disk  D  is  thus  given,  for  z  =  0  and  H  =  0,  by 

E'(M’  u)=  [ £(«:  :!]  =  2  i0’  dkrfAa„{c0kf,  u.  V) 

L  v  J  n.o 


•  C (*,)  •  K a„(k,)  (5) 


where  G(Ar,)  is  given  by  [1 1],  [12] 


£(*,)  = 


-^-(1-/?™)  0 
0 

2 kz 


and  the  reflection  coefficients  for  the  TM  and  TE  polarizations 
are  given  by 

itlkz  +  tklz  tan  (klzd) 


ikz  tan  (kizd)+klz 
ikz  tan 


„  1  /-  d  .  9  \ 

vr  =  -  (  u  —  +  v  — 

h  \  du  dv  / 


and  £z(k„,  u,  v,  z),  u,  v,  z),  £,(*„,  u,  v,  z)  and 

h,  u,  z)  are  the  integrands  of  Ez(u ,  u,  z),  //,(«,  v,  z),  E,(m,  v, 
z)  and  H,(u,  v,  z),  respectively,  ea „(k„)  and  ha,(k„)  are  the 
spectral  amplitudes  of  the  electric  and  magnetic  fields,  respec¬ 
tively,  and  are  obtained  by  equating  the  currents  on  the  disk  to 
the  discontinuity  in  the  tangential  magnetic  field  at  the  disk. 
Thus,  we  get  die  following  relationship: 

*<“■  •>-  [£!£ :!]  -  -s  !>.i*°.<«,*.. «. .) 

L  J  *,a 

■  Ka„{kf),  for  («,  v)  G  D  (3) 


[2 iue  ean(k,)/k„  1 

2 ikt  ha„(kt)/k,  J 


w  W -k\  =  klp- tfVi  -  . 

III.  The  Vector  Integral  Equations  Governing  the 
Current  Distribution  on  the  Disk 

The  current  distribution  on  the  disk  is  governed  by  the 
condition  that  the  current  has  to  vanish  outside  the  disk, 
together  with  the  boundary  condition  on  the  tangential 
component  of  the  eiectric  field  on  the  disk.  These  two 
conditions  give  rise  to  two  dual  vector  integral  equations 
which  can  then  be  solved  using  Galerkin’s  method.  From  the 
condition  on  the  current  we  have 

K(u,  v)  =  0,  for  u>us.  (6) 

From  the  condition  on  the  tangential  electric  field,  we  have 

E,(u,  u)  +  Ef  (u,  v)  =  0,  for  u<u,  (7) 

where  E,(u,  v)  is  given  by  (I),  Ef(u,  v)  is  given  in  Appendix  I 
and  u  =  u,  corresponds  to  the  points  on  the  outer  edge  of  the 
elliptic  disk.  Expressing  these  two  conditions  explicitly  in 
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terms  of  the  current  spectral  amplitudes  we  get 

“S  f  dkpSlan{c0kp,  u,  u) 

n.a 

•  Ka„(A:p)  =  0,  foru>u,  (8) 
2  dkpMa„(c0kp,  u,  v)  ■  G(kp)  ■  Ka,(kp) 

*.a 

=  |  dk'Ma^Cok,,  u,  v)  ■  T a,(kp), 

n.a 

for  u<u,  (9) 

which  are  two  dual  vector  integral  equations  governing  the 
unknown  current  spectral  amplitudes  K a„(kp).  Ta„(Ar„)  is  the 
source  spectral  amplitude  given  in  Appendix  I.  To  solve  these 
equations  using  Galerkin's  approach,  the  current  distribution 
on  the  disk  is  expanded  in  terms  of  the  orthogonal  set  of  the 
TM  and  TE  modes  of  the  magnetic  wall  cavity  since  they  form 
a  complete  set  of  basis  functions. 


we  get 

K0,(M  =  -cSr( k„)  2  £>m.r(o,  0,  C0k „)  •  ac*„m  (14) 

fl.a.ffl 

where 

0,  c0kp)  =  \  ’  du  \  dv 
Jo  Jo 

■  h2  St0r  (Cok„,  u,  v)  •  §a„m(u,  v).  (15) 

Equation  (14)  supplies  a  relationship  between  the  current 
spectral  amplitude  K 0,(kp)  and  the  coefficients  of  the  current 
modal  expansion  aa„m.  The  various  components  of 
0,  cakp)  are  given  by  integrals  /2  and  f3  evaluated  in  Appendix 
II  and  c0,( It,)  is  given  in  [13.  sec.  6], 

The  next  step  is  to  find  the  modal  expansion  coefficients 
a«„m.  This  is  done  by  substituting  the  expansion  of  K0 T(kp) 
given  by  (14)  into  the  second  boundary  condition  given  by  (9), 
multiplying  the  resulting  equation  by  k2Sy'  (u,v)  and  integrat¬ 
ing  over  the  area  of  the  elliptic  disk,  one  finally  gets  the 
following  set  of  equations: 


K(u,  v)  =  £  S a„m(u,  v)  •  aa„„,  u<u,  (10) 


2  7)  '  =  (16) 


where  *a„m  is  a  two-component  vector  whose  components  are 
the  unknown  coefficients  of  expansion,  and  >•)  is 

given  by 


5a»*(M,  v) 


d  t  d  * 

—  <ka„(cQka„m,  u,  v)  —  i>a„(c0ka„m,  u,  u) 

du  dv 

d  ,  3  * 

—  \f/a„(,c0ka„m,  u,  v)  -  —  ta„(c0k a„m,  u,  v) 

dv  du 


an 


where 


d&niCo kctnm,  Hj)  =  0 


(12) 


and 


Ja„(cQka„m,  u,)  =  0. 


(13) 


The  first  column  of  $a„m(u,  v)  represents  the  TM  current 
modes  of  the  magnetic  wall  elliptical  cavity  whereas  the 
second  column  represents  those  of  the  TE  polarization.  Upon 
substituting  the  expansion  in  (10)  into  (6),  we  get 


£  jo" dk't  St0r(cok;,u,  v)  ■  k &(*;) 

=  -  2  $<*<"*(“’ w) '  u<u> 
n,a,m 

=  0,  u>us. 

Applying  the  orthogonality  property  (14)  developed  in  [13], 


where 

y)  =  £  i  dkp  c0r(kp) 

,  Jo 
rJ 

’  co*„)  •  G(kp)  •  Snm,r(a,  0,  c0kp) 

M7)  =  £  f  dkpS'JSJy,  0.  c0k„)  ■  T0r(M 
i.  j  0 

where  the  superscript  t  denotes  transposing  the  matrix.  For 
these  equations  to  be  of  practical  use.  the  modal  series 
expansion  has  to  be  truncated  and  these  equations  can  then  be 
recast  in  the  following  matrix  form: 

A  •  a  =  b.  (17) 

A  is  a  matrix  whose  elements  are  the  matrices  A nm,JS(a,  y), 
whereas  a  and  b  are  vectors  whose  elements  are  the  vectors 
aa„m  and  b„m(a),  respectively.  The  matrix  equation  (17)  can 
thus  be  easily  solved  for  the  unknown  coefficients  a. 

In  the  case  of  a  thin  substrate,  the  structure  approaches  that 
of  the  magnetic  wall  cavity  model  and  in  the  limit  when  d  -»  0, 
the  vector  integral  equations  give  a  solution  similar  to  the 
magnetic  wall  model.  Thus  in  the  case  of  a  thin  substrate 
structure,  the  dominant  modes  are  the  TM  modes.  Further¬ 
more,  if  the  operating  frequency  is  around  w'„m,  which  is  the 
unperturbed  resonant  frequency  of  the  magnetic  wall  cavity 
for  the  TM  polarization  (£-polarization),  we  will  find  that  only 
modes  which  have  a  resonant  frequency  of  in  the 
unperturbed  state,  will  have  a  dominant  contribution  to  the 
currents  on  the  disk  and  all  other  modes  will  be  relatively 
negligible.  Thus,  in  the  case  of  a  thin  substrate  and  for 
operating  frequencies  around  w'„m,  a  single  mode  approxima- 
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lion  can  be  employed  in  which  the  TMa<tm  is  the  only  mode  to 
be  considered. 

IV.  Perturbation  Formula  for  the  Resonant  Frequencies 

The  exact  resonant  frequency  of  the  structure,  acting  as  a 
resonator,  can  be  obtained  by  setting  the  determinant  of  the 
coefficient  of  (17)  to  zero,  to  get 


and  therefore 

L  =  i  (  dz  l  dv  h  £*///„,  evaluated  at  u  =  us. 

J  -4  Jo 

From  (1)  and  (2),  and  by  applying  the  boundary  conditions, 
it  can  be  shown  that 


det  [A]=0. 

As  is  clear  from  this  equation,  the  calculation  of  the 
resonant  frequencies  using  this  equation  is  not  a  simple  task. 
However,  in  the  limit  of  a  thin  substrate,  a  perturbational 
approach  can  be  used  to  calculate  the  resonant  frequencies  in  a 
much  simpler  way.  In  this  limit,  the  elliptic  microstrip 
structure  can  be  viewed  as  a  perturbation  of  an  elliptic 
resonator  with  a  perfectly  magnetic  side  wall.  The  resonant 
frequency  shift  of  this  perturbed  magnetic  wall  cavity  is  given 
by  [10).  [II],  [15) 


A  u  =  u/—u. 


L 

4<Vr), 


(18) 


where 

L  =  -  i  f  j  dS  A  •  (Ef  x  H/) 

.15 

and 


V 

where  E,  and  w,  are  the  electric  field  and  the  resonant 
frequency  of  the  unperturbed  cavity,  H/  and  w/  are  the 
magnetic  field  and  the  resonant  frequency  of  the  perturbed 
cavity  (i.e.,  the  open  cavity),  ( WT),  is  the  unperturbed  time 
average  total  energy  stored  in  the  cavity  and  AS  is  the  surface 
area  of  the  magnetic  wall. 

In  the  unperturbed  case  and  because  the  substrate  thickness 
is  assumed  to  be  small,  the  field  components  are  independent 
of  z.  Thus,  the  only  existing  modes  are  the  TM0*m  modes  for 
which  £;  is  the  only  nonvanishing  electric  field  component. 
Let  us  consider  the  perturbation  of  a  specific  TMMm  mode 
whose  unperturbed  resonant  frequency  of  the  magnetic  wall 
cavity  is  w0„m  given  as  the  solution  of  the  equation 


£/,-  --  2  f  dk>k*  *><•(*<> Ml  +*™> 

€  1  J0 

1  M 

cos  {*lt(z  +  rf)} 
cos  (kud) 

H,t  =  2  J!  dk,k,hyAk,Kl+Rn) 

ry 

sin  {*u(z  +  </)} 


'I'yAcok,,,  u,  u) 


sin  ( kl:d ) 

and,  therefore,  the  u -component  of  Hj  is  given  by 


fyMco*,,  u,  v) 


H/u  =  -h  £  [  £  dk„  ~  hy,(ke)(l  +  f?1 


■) 


cos  {*u(z  +  r/)}  d 
sin  (k[Zd)  dv 


'I'yAcok,,  u,  v) 


-iwt  (  dk,  —  eyr(k,)(\  +  /?™) 

Jo  k * 


cos  {k ,.(z 


cos  (k 


:  +  d)\  d  1 

~~  -  iyAc0k„  m,  v) 

zd)  du  J 


and  hence 


L  =  rA*Ja„(Cok,,  us )  S  ^  f0  dk"  h<*Ak» 
■  (\+R™)Jar(c0k,,  us)QdrAcok,,Cok,) 

+  «(’  dk  ley^Kl+Z?™) 

Jo  k, 

tan  ( kud ) 


k\z 


J'ar(C0k,,  u$)Na*Acok,,  c0k,)  (20) 


kotnm,  us)  -  0 

where 

For  this  mode,  the  unperturbed  electric  field  and  resonant 
frequency  are  given  by 

E,-i  A  4>a„(Cok,,  u,  v),  k,  =  ka„m  (19a) 

and 


where  Qdm(c0k,,  c0k,)  and  Sanr(cok„  c0k,)  are  defined  and 
given  in  Appendix  II  and  by  a  we  mean  the  parity  opposite  to 
a. 

In  the  limit  when  d/a  —  0,  and  from  (10),  (19),  (4),  (14), 
and  the  results  of  Appendix  II,  aanm,  har(k,)  and  ea,(k„)  can 
be  approximated  as 


,  *v4  w.e,  1 

=  —  k„  car(k,)  —  r 

2  <**  k;-k] 

J&h (c0k, ,  us)J  <*Ac0kp,  u j ) Afat „r(c0k, ,  Cgk, )  (21) 


( 


®  Chorum 


(19b) 
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hdr(k„)  =  —-  cd,(k„)  -1 - 

2  k.  ui,no 

■  Jd,(c0kp,  u,)Joc„(Cok,,  us)Qct„(c0kn  c0kp).  (22) 

Finally,  the  unperturbed  time  average  total  energy  stored  in 
the  cavity  is  given  by 

«il'4|J  d  f  dv  (  1  du  h:[ict„(c0k,,  u,  i>)]2 
2  Jo  Jo 

and  from  Appendix  II,  we  obtain 

1  \A  I2 

( W'V),  =  --  re,d  ——  Ja„(c0k,,  u,) 

4  k, 

■  I  J'an(c0kp,  us) I  Ra„(c0k,).  (23) 

Substituting  (20),  (21),  (22),  and  (23)  into  (18)  we  finally 
get  the  required  perturbational  formula  for  the  resonant 
frequency  shift  of  the  TMa„m  mode. 

V.  Calculation  of  the  Input  Impedance 

The  input  impedance  across  the  terminals  of  the  probe  P  is 
given  by  [12),  (16) 

Zm=  jd^E,  •  J 


where  E(  is  the  total  electric  field  in  the  substrate  region  due  to 
the  current  on  the  disk  D  as  well  as  the  current  on  the  probe  P. 
The  volume  integration  in  the  expression  of  the  input 
impedance  is  carried  out  over  the  volume  of  the  probe.  J  is  the 
current  distribution  on  the  probe.  To  simplify  the  analysis,  the 
current  distribution  on  the  probe  is  modeled  as  a  uniform 
cylindrical  current  sheet  of  radius  a  and  hence  J  is  given  by 


thus,  the  input  impedance  is  given  by  [12) 


z«=r~  S’ dk.  r  y«(M) 

f  2  ( 2k\e*izdn  sin  (kl;d/2) 
■  (I  -  /?™)  tan  (kud/ 2)  e"kizd/i 


)  k1 

T‘t 


/  6[\p-p0\-a\ 

j  =  i - - — 

2r  |p  — Pol 


-d<z< 0 


where  po  is  the  vector  connecting  the  origin  of  the  coordinate 
system  to  the  center  of  the  probe.  The  input  impedance  is 
hence  given  by 

-Ji  {  dV(E’t  +  E\z)Jt 

where  Eu  is  the  electric  field  due  to  the  current  on  the  probe. 
E\t  is  the  electric  field  due  to  the  currents  induced  on  the  disk 
and  can  be  shown  to  be  given  by 

E‘u~  ~ 2iut,  ^  L 


fl  +  2  — +— i— 
l  X™  JJ  2iue,I 


2  f  Jo(kpa)[Ky pik,,)],, 

Z  0  k 

(1  +  /?™)  tan  (kltd)+yp(c0k„  u0,  u0) 


„  TM_  fl  kz 

A  in  — - 


[Ky^Ar,)),,  is  the  u-component  of  Kyp(k„)  and  is  given  by  (14) 
and  (16).  In  the  limit  of  a  thin  substrate,  and  since  the 
dominant  modes  are  those  of  the  TM  polarization,  the  u- 
component  of  the  current  on  the  disk  can  thus  be  shown,  from 
(14),  (15)  and  the  results  of  Appendix  II,  to  be  given 
approximately  by 

[Ky /»(*„)]„  =  r  cyp(k,)J'yp(c0kp,  us) 

T-l  TM  (toffs)1  .  .t  . 

'S  kl-ik^)*™01  yr”  ^ 

■  Nyrpicoky c0k, )  (24) 

where  ay™  are  the  TM  current  expansion  coefficients.  Under 
the  single  mode  approximation  (corresponding  to  the  TMa„m 
mode),  we  thus  obtain 

[Kctp(kp)]u  =  T  cap(k0)J'ap(cok„  u ,) 

tm  (.ka„my  e 

■  aapm  —  ■  —7 — —  Ja„(c0kanm,  u,) 

*  p  \K&nm) 

N<Xnp{Cokcinm ,  c0kp)  (25) 

and  hence  the  input  impedance  is  given  approximately  by 


r  dk>  it  yo(**fl) 

4to*.  Jo  k,. 


[it 


2k\eik'tdn  sin  (kHd/ 2) 


cos  {<:itU-Fd)} 
cos  (kud) 


(1  +  R™)*yP(Cok„  u,  v) 


+  i  rf  (1  -  A?™)  tan  (kud/ 2)ettnrf/2 

*it 


1 


7 

! 
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■  £  cctp(kp)Nctnp(c0kct  n  m  *  c0kp) 


TM  .  k‘ 

aa,„  Ja„(c0ka„m,  us)  dkp  — 2  J0(kpa)(  1  +  RJH) 

J0  A|j 

(kapm)  ^  ca^k^J'apicok,,,  us) 


'  J'otp(.c0kp ,  us)Sctp(c0kp,  v)Ha^,(c0kpt  u). 

In  the  far  field,  it  can  be  shown  that  c0  cosh  (u)  -*  p  and  that 
the  asymptotic  expansion  of  Ha^(c0kp,  u)  is  thus  given  by 


tan  (kud) 


k2p-(ka„m)2 

•  Nanp(c0ka„m,  c0kp)\l/ap(c0kp,  «o.  u0). 

VI.  Radiation  Pattern 

In  this  section,  the  far-field  expression  for  the  electric  field 
in  the  thin  substrate  limit  will  be  developed.  In  this  case  and 
because  the  structure  approaches  the  magnetic  wall  cavity 
model  (where  the  probe  becomes  almost  shielded  from  the 
outside  region),  the  contribution  of  the  current  on  the  probe  to 
the  radiation  field  is  negligible  compared  to  that  of  the  currents 
induced  on  the  disk. 

From  (I)  and  (4),  the  z -component  of  the  electric  field  due  and 
to  the  current  on  the  disk  in  the  upper  half-space  is  given  by 


Ha{''(c0kp,  u)~  Lctp{Cok,,) 


1 


sjrkTp 

■  exp 


\)f,Lap(c0kp)H{0"(kpp) 


where 


Lep(c0kp)  =  Sep(c0kp<  0) 


Lop(c0kp)  =  S'Op(e0kp,  0). 


Thus,  in  the  far  field,  i.e..  for  large  values  of  p  and  z,  the 
integrand  is  rapidly  oscillating  and  the  integral  can  hence  be 
TM  t  evaluated  using  the  method  of  stationary  phase.  A  stationary 
w7p(c0A„,  u,  v)(l  -  R  )e'  exjsts  at  kp  =  k  sin  (0).  where  0  =  tan-1  (p/z).  Thus, 


Under  the  single  mode  approximation  (for  the  TMa„m  we  8ct 
mode),  [Kap(Ar„)]u  is  given  by  (25)  and  hence  £.  is  given  by 


”  rr  QQnm  kotnmt  wr)  i  dkpk ^ 

2  lux  Jo  p 


r, - —  [l 

4  iux  -y  2 


( ka „„)2 


(1  -R™)e'*t: 


k2  sin  (0)  cos  (0) 

(ka„m)2 


k\-(ka„„)2 

■  Y,  Cotpik^Na^Coka  nm »  Cokp) 

P 

■  J'apicokp,  us)Sctp(c0kp,  v)Jap(c0kf,  u ). 

It  can  be  shown  that  the  characteristic  value  b„  together  with 
the  expansion  coefficients  D"m(c0kp)  and  F"m(c0kp)  of  the 
Mathieu  functions,  defined  in  (13,  appendix  I],  are  even 
functions  of  k„.  Thus 

J'atp(-c0kp,  u,)  =  (-l)i>J'ap(Cokp,  us ) 

Ha"'(Cok,e‘\  u)=  -  (-  \)»Ha™(c0kp,  u) 

and  Na„m(c0ka.m,  c0kp),  Sap(,c0kp,  v)  are  even  in  kp.  whereas 
cap(k,)  is  odd  in  kp.  Hence,  the  integral  along  the  positive  real 
axis  in  the  expression  of  £,  can  be  extended  to  an  integral  over 
the  whole  real  axis,  to  get 

x  t 

tz  =  4hx  aa"m  Ja^coka^<  u,) 

■  \\dkpki7!k°::]\j\-R^)e^2 


^8f/i(Co^Otfiffi,  Mj) 


(ka„m)2~k2  sin2  (0) 

■  [1  -R™(kp  =  k  sin  0)]  £  (-  1  y  cap{k  sin  0) 

p 

■  Lap(c0k  sin  6)Na„p(c0ka„„,  c0k  sin  0) 
J'ap(c0k  sin  0,  us)Sap(c0k  sin  0,  <i>) 

■  dkp^e'^'H^kpp). 

J  -  «  A- 

Using  Sommerfeld's  identity  (14),  we  finally  obtain 


—  E 

:  2  yj  2 


aa™k2  sin  (0)  cos  (0) 
(ka„m)2 


( ka„m)2-k 2  sin2  (0) 


Ja„(c0ka  nm  *  ^i) 


9tkr 


k\-{k^m)2 


[\ -  R™(kp  =  k  s\n  0)]  —  £  (-DpcaP(*  si 

r  P 

Lap(c0k  sin  6)Na„p(.c0kanm,  c0k  sin  0) 

•  J'ap(c0k  sin  0,  u,)Sap(c0k  sin  0,  <t>) 


sin  i 


'  \ 


* 
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where 

r  =  \lpl  +  z2 

and  in  the  far  field  v  ~  <t>.  Similarly,  the  expression  for  Ht  in 
the  far  field  is  given  by 

Ht - -  aa™  k  sin  (6)Jan(c0kanm,  us) 

■  [l  +  £TE(*,  =  *sinfl)l  —  •  V  (-Wcd^k  sin  8) 

r  ** 

p 

■  Ldp(c0k  sin  8)Qa„p(c0ka„m,  c0k  sin  8) 

•  Jd„(c0k  sin  8,  u,)Sdp(c0k  sin  8,  4>) 

where  d  =  e  (even)  if  a  =  o  (odd)  and  a  =  o  (odd)  if  a  =  e 
(even).  Finally,  the  8-  and  ^-components  of  the  electric  field 
are  given  by 

£«  =  -£*/sin  (8),  and  £,=  -  v^,/#  //./sin  (8). 


and  g(r,  r')  is  the  solution  of  the  following  differential 
equation 

(VJ  +  *J)g(r,  r ')=  &(u-u')6(v-u')6(z-z'). 

From  the  results  of  [13),  g( r,  r')  can  be  expanded  in  terms 
of  the  scalar  elliptic  harmonics  as  follows: 

g(T,r')=‘-Y.  (  dk„  wa„(kf)  ^(Cok,,  u’,  v') 

2  “  0  k* 

■  ia„(c0kp,  u,  uje'**1*-  '1 
and  hence  the  electric  field  is  given  by 

E(w.  i>,  z)=--^2(I*l  +  W)  (  dr' i(r')  ("  dk, 

7;  Jvp 

•  ~  wan(kf)ian(c0kf.  u' ,  u,  v)e,ktu't"[ 

k. 


VII.  Conclusion 

The  scalar  and  vector  Mathieu  transforms,  developed  in 
[13],  allow  the  rigorous  formulation  of  the  elliptic  disk 
microstrip  antenna  which  has  long  been  thought  almost 
impossible  to  formulate  exactly  [6]. 

It  is  shown  that  the  current  distribution  on  the  disk  is 
rigorously  derived  using  these  transforms  and  that  it  is 
governed  by  vector  integral  equations.  These  vecror  integral 
equations  are  then  solved  using  Galerkin's  method  in  which 
the  current  distribution  on  the  disk  is  expanded  in  terms  of  the 
TM  and  TE  current  modes  of  the  magnetic  wall  cavity  model, 
which  form  a  complete  set  of  basis  functions.  The  limit  of 
small  substrate  thickness  is  then  applied  to  these  equations  to 
get  simple  approximate  expressions  for  the  current  ampli¬ 
tudes. 

The  resonance  in  the  elliptic  disk  structure  is  then  analyzed 
using  two  different  approaches:  Galerkin's  method  and  a 
perturbative  approach.  The  input  impedance  together  with  the 
radiation  pattern  are  derived  both  exactly  and  in  the  small 
substrate  thickness  limit. 

Appendix  I 

Field  Expressions  Due  to  Probe  Excitation 

First,  we  will  develop  field  expressions  due  to  the  probe 
excitation  in  the  unbounded  medium.  In  the  next  part  of  the 
appendix,  we  will  then  generalize  the  expressions  to  include 
the  effect  of  stratification.  The  electric  field  E  is  governed  by 
the  following  differential  equation 


where  Vp  is  the  volume  of  the  probe.  Thus, 

•  wa„(kp Wa* (cakp,  u0,  uQ))^a„(c0k„,  u,  \>)f(k„  z) 
where  /(*„,  z)  is  given  by  [12] 

/(*,,  ?) 


'gdglt-t'l 


f 


cos  {k^z  +  d/2)}  -  1], 

itcl 

2 

eik'U*d/i\  sin  (ktd/2\ 

"t 


-d<z<~0 
0 <z<  -d 


Thus,  in  the  case  of  the  structure  considered  in  Fig.  1  and  by 
applying  the  stratified  medium  formulation  [14],  the  z- 
component  of  the  electric  field  in  the  air  region  and  the 
substrate  is  given  as  follows  [12]. 

In  the  air  region  (z  >  0): 


£f(«.  v.  z)=  ^  (  dk>Ti  Jo(kfa)wan(kp) 

•  (1  - R™)4ia„(c0kp,  u0,  v0)ia„(c0kp,  u,  v)eUctz.  (26a) 


VxVxE(u,  u,  z)-k2E(u,  u,  z)  =  iwnol(u,  v,  z)  In  the  substrate  (-d  <  z  <  0): 


where  J(u,  v,  z)  is  the  current  distribution  on  the  probe.  We 
will  represent  the  solution  of  this  equation  in  terms  of  the 
dyadic  Green's  function  in  the  elliptic  coordinate  system, 
which  has  the  following  form 

G(r,  r')=  (l  +  p  $(*.  *') 


£*>,  v,  z) 

=—  y  \  dk„-~-  J0(k,a)waAkt) 

<*’  Z  Jo  ku 

■  ta„(c0kp,  w0.  voWa,(Co*e,  «.  ”) 
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k2e'ktidn  cos  {kXt(z  +  d/2)}  -k\ 


ik\  sin  (k,zd/l)e'kud/2 

(1 


[{l+R\?e‘ku*) 


■d/2). 


(26b) 


where  /?™  and  are  given  in  Sections  II  and  V, 
respectively.  At  z  -  0,  the  transverse  component  of  the 
electric  field  is  given  by 


E f(u,  o)  = 


£•;(«, o) 
£>•  v) 


} 


=  J  dk,Man(c0kQ,  u,  o)  ■  T a„(*„)  (27) 


where 


and 


Ta„(*,)=  [Pafk)] 


Pa,,(k,)  =  ^-  -f  k0Mkta)waH(kf) 

2“* 

•  (1  -/?™)^a,,(Co*„  M0.  wo)- 
Appendix  II 

Evaluation  of  /,  =  j J’  dv '  j;;  du'Zi2  V'oufco*,,  m',  o') 
■  WrlCok u',  o'):  +ct„(c0ke,  u',  o')  and  t0r(cok’,  m',  o') 


are  governed  by  the  following: 

du7^  ^a"  +  ^a" +  *2*2  ^a" =  0  (28) 

^Mr+^Mr  +  h2k'2Wr  =  0.  (29) 

du  2  do  2 

Multiplying  (28)  by  ^0,  and  (29)  by  ^a„  and  subtracting, 
we  get 

d  d 

V'Cr,  —  \l>0r-\k0r  —  \pa« 

dM  du 

d  d 

Wr-Wr  — 
do  du 

+  **(*,'  2  -  * J)  ^O,  i0r  =  O. 

Integrating  over  m'  from  0  to  u  and  over  o'  from  2  to  2x. 
Using  the  periodicity  property  of  the  angular  functions 
together  with  the  fact  that  J'ym(c0kp,  0)  =  0  and  finally 
employing  the  orthogonality  relationship  of  the  sine  and  cosine 
functions,  we  get 

7I  =  rr_T-2  [Jot*(cok,,  u)J’arlc0k’,  u ) 

P  P 

-Ja„(Cokp,  u)Jar(c0k^  ,  «)]iVa„(c0<:,,  c0k^  )6^ 

provided  that  k„  *  k't  and  both  n  and  r  have  to  be  even  or  odd 
(/i  vanishes  if  n  and  rare  of  different  parity).  Samr(c0kc,  c0k p 
is  defined  as  follows: 

1  r*» 

Na„{c0k„,  c0Ar  ')  =  -  l  dv  Sa„(c0k,,,  v)Sar(c0k' ,  o) 

f  J0 
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hence, 

Ne„e(c0k„  c0kp  =  J)'  (1  +  6m0)D’,m(e0kll)D'm(cok 
No„(cak„  c0k' )  =  ]£'  FHm(c0kl,)FrJcok'). 

m 

If  k'  -*  Ihe  following  two  cases  have  to  be  considered. 

Case  (i):  r  *  n:  In  this  case  and  as  k'a  —  kp,  the  quantity 
between  the  square  brackets  in  the  expression  for  /,  has  a 
nonvanishing  value  whereas  Na „r(Cokp,  c0k^)  vanishes  (see 
[13,  eq.  (4)]).  Applying  L’Hopital's  rule,  we  obtain 

x 

h  ~  [/a,(c0Arfl,  u)J'ar(c0kp,  u) 
-J'ap(c0kp,  u)Ja,(c0k ,,  «)] 


NctnrfCok, 


Case  (ii):  r  =  n:  In  this  case  and  in  the  limit  when  k‘t  -* 
kp,  Napp(c0kp.  c0k'p)  -  Ra„(.c0k„)  (see  [13,  eq.  (3)1)  whereas 
the  quantity  between  the  square  brackets  vanishes  and  hence, 
we  get 

[s;  ■**.<**;. 

Rct„(c0  kp  . 

Evaluation  of  I2  =  du'  jj'  dv'[d/du'  i£a„(cokfl,  u',  t/') 
d/du'  il/0r(cok^  u',  u')  -  d/dv '  \(ia„(c0kp,  w') 

d/du‘  MriCoky  u' ,  u')J.  Using  integration  by  parts,  two 
alternative  forms  for  /2  can  be  obtained.  The  first  is 


r2  =  *J'ap(c0kp,  u)Jar(c0k'  u) 


■  N<x„(c*kt,Cok'f)ba6  +  k\Ix 


and  the  second  is 


/2  =  ■wJa„(Cokt,  u)J'<x,(c0k'e ,  u) 

■  Na„(c0kp,  Cok'^  +  k'2^. 

n  and  r  have  to  be  of  the  same  parity  (/2  vanishes  if  they  are  of 
opposite  parity). 

Evaluation  of  /j  =  jj  tfu'  (J*  dv'[d/du’  \(ia„(Cokp,  u',  o') 
d/dv'  t0,(cok',  i/')  +  d/dv’+an(c0ke,  u',  v') 

d/du' 4i0f(cok  't,  u' ,  o')].  Using  integration  by  parts,  it  can  be 
shown  that  /}  is  given  by 

fjs  —  e  Jotp(.Cok0,  u )  Jdr(c0 k ^  ,  u)Qaw(to^i  Co k  ^  )<5 ^ 

where  3  =  e(even),  if/3  =  o (odd)  and /3  =  o(odd),  if/3  =  e 
(even),  and 

1  rJ* 

c0k  )  =  -  I  duS  a„{c0kp ,  u)Sdp(cok  ,  v) 

5T  J0  * 


QeeAcakp,  c0k;)=  mD”Jcokp)F'Jc0k;) 

m 

Qo„,(c0kp,  c0k;)='£i’  rnF"m(c0kp)D'm(c0k'o) 


with  n  and  r  of  the  same  parity  (i.e.,  both  n  and  r  have  to  be 
either  even  or  odd). 
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Transient  Response  of  a  Vertical  Electric  Dipole  (VED) 
On  a  Two-Layer  Medium 


S.  Y.  Poh  and  J.  A.  Kong 
Research  Laboratory  of  Electronics 

Department  of  Electrical  Engineering  and  Computer  Science 
Massachusetts  Institute  of  Technology 
Cambridge,  MA  02139,  USA 

Abstract-  The  transient  electromagnetic  radiation  by  a  vertical  electric  dipole  on  a  two- 
layer  medium  is  analysed  using  the  double  deformation  technique,  which  is  a  modal  tech¬ 
nique  based  on  identification  of  singularities  in  the  complex  frequency  and  wavenumber 
planes.  Previous  application  of  the  double  deformation  technique  to  the  solution  of  this 
problem  is  incomplete  in  the  early  time  response.  In  this  paper  we  show  that  the  existence 
of  a  pole  locus  on  the  negative  imaginary  frequency  axis,  which  dominates  the  early  time 
response,  proves  crucial  in  obtaining  the  solution  for  all  times.  A  variety  of  combinations 
of  parameters  are  used  to  illustrate  the  double  deformation  technique,  and  results  will  be 
compared  with  those  obtained  via  explicit  inversion,  and  a  single  deformation  method. 


I.  INTRODUCTION 

The  time  response  of  a  VED  over  a  two-layer  nondispersive  dielectric  was  ana¬ 
lyzed  by  Ezzeddine  et  al.  [1]  employing  the  geometrical  optics  approach  for  early 
arrivals,  together  with  an  explicit  inversion  scheme,  analogous  to  the  C&gni&rd-de 
Hoop  method  [2,3],  that  is  valid  for  all  times.  The  solution  is  expressible  as  single 
integrals.  As  the  latter  technique  is  inapplicable  for  problems  involving  dispersive 
media,  Ezzeddine  et  al.  [4]  further  analyzed  the  pulse  response  of  a  VED  over 
a  two-layer  medium  by  applying  the  double  deformation  technique  [5]  which  in¬ 
volves  complex  plane  deformation  in  both  the  wavenumber  and  frequency  planes. 
The  resulting  solution,  although  numerical  in  nature,  does  not  require  excessive 
computation,  and  may  be  readily  extended  to  consider  dispersive  media.  The 
results  as  presented  in  [4],  however,  lacked  completeness  insofar  as  the  early  time 
solutions  for  both  lossless  and  dissipative  media  are  concerned. 

In  this  paper,  we  seek  to  resolve  the  problem  encountered  by  Ezzeddine  et  al. 
[4]  in  the  early  time  response  for  the  radiation  of  a  VED  on  a  two-layer  medium. 
The  existence  of  an  additional  pole  locus  on  the  negative  imaginary  frequency 
axis  proves  crucial  in  improving  the  early  time  response.  The  analysis  presented 
is  valid  for  both  non-dispersive  and  dispersive  media.  However,  the  emphasis  in 
the  presented  results  is  on  non-dispersive  media  since  much  of  the  effort  and  moti¬ 
vation  is  to  complete  an  initial  picture  regarding  the  double  deformation  approach. 
Some  variety  in  the  choice  of  parameters  is  used  to  illustrate  the  double  defor- 
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mation  technique,  and  results  will  be  compared  with  those  obtained  via  explicit 
inversion,  and  a  single,  wavenumber  plane,  deformation  approach. 


n.  FORMULATION 

The  geometrical  configuration  to  be  considered  [Fig.  1]  comprises  a  vertical  electric 
dipole  (VED)  located  at  the  origin  on  the  surface  of  a  two-layer  medium  and  the 
observation  point  a  distance  p  away,  on  the  same  horizontal  plane.  The  dipole  is 
assumed  to  be  excited  by  a  current  pulse  which  is  zero  for  t  <  0  and  a  damped 
sinusoid  for  t  >  0 ,  e.g. 

t(<)  =  lo  tn  sinu»o<  e-Ol0<  u_i(t)  for  n  =  0, 1,2, ...  (1) 

where  u_j(<)  is  the  unit  step  function.  This  choice  of  function  allows  us  to  simu¬ 
late  a  variety  of  waveforms  by  varying  the  values  of  u>o  and  ao-  The  corresponding 
Fourier  transform  is  given  by 


where  k0  =  u i/c  and  c  is  the  speed  of  light  in  free  space. 

The  upper  half-space  medium  is  assumed  to  be  free  space  with  permittivity 
e0  and  permeability  p0-  Using  cylindrical  coordinates  p,  <j>,  and  z,  the  time- 
harmonic  response  for  the  magnetic  field  at  a  point,  (p,$,z),  in  the  upper  half 
space  is  given  by  [6] 


where 


H  =  = 


Sir  ,/SIP 
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dk ,  V) 


*01 


1  + 


s™ 


rTM  *01  +  Rn«l2kud 
1  + 

D  _  Srlfeo*  ~ 

*  &rlkoz  +  *lz 

n  _  ~  Srl*2z 

12  ir2*l*  +  £rl*2z 


(3) 

(4) 


klx  +  =  *o  = 

*?,  +  kp  =  =  ^Poirl^o  —  ir\k0 

k2i  =  -2  =  u,2Moir2<o  =  ir2 ko 

and  SIP  is  the  Sommerfeld  integration  path.  For  the  configuration  of  interest  [Fig. 
1],  z  =  0. 
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Figure  1.  Geometrical  configuration  of  the  problem. 

In  the  time  domain,  the  magnetic  field  may  be  expressed  as  the  inverse  Fourier 
transform 

H^t)  =  ifle|^°°dWe-"<^(u;)| 

The  evaluation  of  the  resulting  double  integral  is  difficult  both  analytically  and 
numerically  [7-8].  In  the  double  deformation  method  [5],  the  paths  of  integration 
are  deformed  to  the  steepest  descent  paths  (SDPs)  in  both  the  kfi  and  u>  planes. 
The  double  integral  corresponding  to  the  transient  magnetic  field  is  then  reduced 
to  single  integrals  associated  with  the  contributions  of  singularities  enclosed  in 
the  path  deformations  and  double  integrals  along  the  steepest  descent  paths  in 
the  kp  and  w  planes. 

m.  DEFORMATION  IN  THE  k p  PLANE 

For  the  integrand  in  the  kp  plane,  we  find  two  branch  points  at  kp  =  k0  and 
bp  =  k_2-  There  is  no  branch  point  at  kp  =  &i»  corresponding  to  k\t  =  0,  since 
the  integrand  in  (3)  is  even  in  k\x.  We  first  deform  the  original,  Sommerfeld 
path  of  integration  (SIP)  in  the  kp  plane  to  the  vertical  steepest  descent  paths, 
SDPl  and  SDP2  [Fig.  2]  about  the  branch  points  at  k0  and  £2  respectively.  The 
branch  cuts  as  shown  in  Fig.  2  are  chosen  such  that  kox  and  kix  have  positive 
real  parts  on  their  respective  upper  Riemann  sheets,  on  which  SIP  is  defined,  and 
negative  real  parts  on  their  lower  sheets.  In  the  process  of  deformation,  we  find 
that  the  deformed  path  traverses  three,  namely  the  sheets  corresponding  to  upper 
sheet  for  both  kox  and  Jkj*  (UU),  lower  sheet  for  kox  and  upper  for  &2z  {LU) 
and  lower  sheet  for  both  kox  and  *2*  (££),  of  the  four  possible  sheets.  In  our 
analysis,  we  choose  to  fix  the  real  part  of  kox  and  of  k^x  as  positive  regardless  of 
Riemann  sheet  but  apply  the  appropriate  signs  explicitly  in  the  integrand  instead. 
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In  other  words,  from  (4),  on  the  Riemann  sheet 

dTM  _  *01  +  R\2*'2kud  dTM 
Ruu  ~  1  + 


denoted  by  {UU), 

Re(koz),  Re(k2z)  >  0 


while  on  the  sheet  ( LU ), 


R 


TM 

LU 


1  -  R<nR\2el2kud  =  1/aTM 
Roi  -  Ri2«l2kud 


Re(koz),  Re(l b2x)  >0 


since 


=  1/*01 


(5) 

(6) 


On  the  sheet  ( LL ), 


fl.i2  -  R0ie'2k"d 
*01*12  f  C‘**l*^ 


Re(fcot),  /?e(Jb2r)  >0  (7) 


because 

*12!.  ,  =  1  /  *12 
*2i  =  -*2x 

In  the  A.-p  plane,  pole  singularities  may  also  be  present,  on  the  traversed  Rie¬ 
mann  sheets,  as  zeros  of  the  denominators  of  the  reflection  coefficients  in  (5) -(7). 
The  positions  of  these  poles  vary  as  a  function  of  frequency  u/. 


Figure  2.  Original  and  deformed  integration  paths  in  kp  plane  with  en¬ 
closed  poles. 

o  :  pole  intercepted  on  the  sheet  Re(koz)  >  0;  Re(k2z)  >  0. 
x  :  pole  intercepted  on  the  sheet  Re(koz)  <0,  Rt{k2z)  >0. 

No  poles  intercepted  on  the  sheet  Re{koz  )  <  0,  Rt{k2z)  <  0. 

At  this  point,  the  original  kp  integral  is  replaced  by  contributions  from  any 
intercepted  singularities,  the  integration  over  the  semi-circular  arc  and  the  inte¬ 
grations  along  the  vertical  paths  SDPl  and  SDP2. 
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By  Jordan’s  lemma,  the  integration  over  the  semi-circular  path  evaluates  to 
zero.  We  note  that  along  SDP1,  kp  —  k0  +  iq ,  where  q  >  0  and  the  contribution 
to  the  time-domain  magnetic  field  is  given  by 


=  ^4 R*  (  r  dk0e-ik*TI{ Jb0)  [“dq-J^L 

[JO  Vo  qy/kl-(k0^ 


H[l)  ((fc0  +  iq)p)  [2  +  ^ 


and  along  SDP2,  kp  =  k2  +  *?»  9  >  0, 

BtsDPi(r)  =  ^R'  |  jT ik"-**rl(l")  f‘ 


{Jc2  +  iq )2 

'kl  -  (k2  +  iq P 


^  ((^2  +  *«)P)  *,-*,+»*  } 


where  we  have  performed  the  normalization  k0  =  m/e,  t  —  ci  and  I(u)  = 
I0I{k0)/c,  and  where  the  square  roots  are  defined  as  having  positive  real  parts. 

Therefore  after  deformation  in  the  kp  plane  the  total  magnetic  field  H^(t)  is 
given  by 

S4>(r)  =  HjSDPlir)  +  B4SDP2(T)  +  ^  B+MUu(T) 

M 

+  5Z^Afi^(r) +'52h*mLL(t)  (io) 

M  M 

where  H^sDPl(T)  “d  ^5DP2(T)  8re  the  contributions  from  the  steepest  de¬ 
scent  path  integrals  and  ff^\fuu(T)y  H*MLu(T)  aQd  are  the  contri¬ 

butions  from  the  intercepted  poles  in  the  deformation  through  the  three  Riemann 
surfaces. 

To  obtain  the  contributions  from  the  intercepted  poles,  we  have  to  first  deter¬ 
mine  their  positions  in  the  complex  kp  plane  as  a  function  of  frequency  ka.  We 
first  write 

*TM  .  flVaM  (iu) 

*TM  ,  (Ui) 

w  Svw 

rW  =  mo 


where 


LL  rf*  (kp,k0) 

9t^(kp,k0)  =  1  +  J2<ji  J?i2e’2*u<* 
f?{kP,ko)  =  +  Ri2ei2kud 


» 
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r,M(kp,*o)  =  *01*12 +  «‘2W 
=  *12  +  *01«*2*1*^ 

To  determine  the  existence  and  positions  of  the  poles,  we  have  to  solve  for  the 
zeros  of  g^(kp,ko),  f?{kp,kc)  and  r^(fcp,  k0)  in  the  complex  kp  plane.  The 
residues  of  these  singularities  contribute  to  the  field  solution  only  if  they  lie  in 
the  respective  regions  of  interest.  For  instance,  for  the  zeros  of  9^{kptk0)  to 
contribute,  they  have  to  lie  in  the  ( UU )  region  of  the  kp  plane  [Fig.  2]. 

Denoting  the  zeros  of  g^(kp,k0)  by  K\(uu,  the  zeros  of  /,M(fcp,  k0)  by 
KmLU>  and  those  of  r^(kfi,k0)  by  K mil,  we  have  that  kp  =  KmuV  satisfies 

ln(J2oi*i2)  +  i2k\td  —  i(2M  +  l)ir  Af  =  integer  (12) 

while  kp  =  Kffiu  satisfies 

In  j  +  i2kud  =  i{2M  +  l)ir  AT  =  integer  (13) 

and  kp  =  K M LL  satisfies 

—  ln(f?oilii2)  +  i2k\zd  =  i(2M  +  l)ir  Af  =  integer  (14) 

where  all  Re(y/o)  >  0.  The  pole  loci  for  K^uu  and  Kmhj,  normalised  to  layer 
thickness  d ,  for  varying  Af  are  shown  in  Fig.  3  for  =  3.2  (ice),  and  e2  =  80 
(water).  Details  on  the  determination  of  these  loci  are  presented  in  [9].  Although 
there  axe  two  sets  of  poles  for  each  Af  owing  to  the  modal  equations  being  even 
in  kp ,  we  have  shown  only  the  set  that  lies  in  our  domain  of  interest.  As  noted 
in  Fig.  2,  we  find  no  solution  for  Kmh  corresponding  to  (14)  in  the  region  of 
interest  which  is  designated  ( LL ).  This  will  be  the  case  whenever  «2  >  el»  and 
is  explained  later  in  this  section. 

Both  K\juu  and  Kj^lu  converge  towards  kp  =  k\  for  kad  — ►  oo  but  kp  = 
fcj  is  not  a  pole  singularity  of  the  reflection  coefficients.  Mathematically  this  is 
due  to  the  fact  that  R^if  and  *jT^  remain  finite  for  this  value  of  kp. 

From  the  deformation  process  illustrated  in  Fig.  2,  we  find  that  the  singularities 
of  importance  to  the  left  of  the  Re(Jbp)  =  line  are  those  {Kmuu)  due  to 
R$j-  since  the  path  remains  on  the  top  Riemann  sheet  intercepting  only  the 
corresponding  poles.  On  the  other  hand,  to  the  right  of  the  Re{kp)  =  k0  line  and 
to  the  left  of  the  *e(fcp)  =  Re(ki)  line,  the  deformed  path  crosses  the  branch  cut 
into  the  lower  Riemann  sheet  of  k0t  but  remains  on  the  top  Riemann  sheet  for 

mw 

k2z  and  therefore  intercepts  poles  {Kmlu)  attributed  to  Rfa  .  Furthermore 
we  find  that,  for  varying  1 fe0,  the  pole  loci  for  Kmu  intersects  the  Re(kp)  =  k0 
line  at  two  points  corresponding  to  k0  =  and  k0  =  k^ 2-  This  explains  why 
we  need  only  integrate  over  the  range  k&g  j  to  fcj/j  for  k0  when  we  consider 
the  field  contribution  of  Kmjju.  For  the  Kmlu  l°ci»  *be  intersection  with  the 
Re(kp)  =  k0  line  occurs  only  at  one  frequency  corresponding  to  k0  =  and 

there  is  no  intersection  with  the  Re(kp)  =  Reiki)  ^ne- 

Physically,  A 3^3  corresponds  to  the  cutoff  frequency  of  the  Af  th  transverse 
magnetic  (TM)  mode  guided  in  the  bounded  layer.  The  integration  from  k^ 3  to 
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Figure  3.  Pole  loci  oa  kp/k0  plane  for  =  3.2,  «2  =  80,  and  positive  k0. 

(a)  Pole  loci  (Kfguu),  M  =  1,2,  3.  Region  of  interest  is  left 
of  Re(kp)  =  k0.  (b)  Pole  loci  ( Kmlu )*  M  =  1,  2,  3.  Region 
of  interest  is  right  of  Re(kp)  =  k0  and  left  of  Re(kp)  =  Re(k2) 
(not  shown). 

oo  sums  the  contribution  of  the  Af  th  guided  mode  over  all  frequencies  of  prop* 
agation.  We  End  that  the  excited  Kauu  poles  correspond  to  leaky  wave  modes 
in  both  the  upper  and  lower  half-spaces  having  k0  >  Re(kp)  >  0,  Im(kp)  > 
0;  Re(kox)  >  0,  Im(kox)  <  0,  and  Re(k2X)  >  0,  Im(k2x)  <  0,  while  the  excited 
KMLU  poles  correspond  to  turface  wave  modes  having  Re(kox)  <  0  and  Im(k0t ) 
>  0,  in  the  upper  half-space  and  leaky  waves  in  the  lower  half-space.  We  comment 
that  any  K^n  pole  singularities  would  correspond  to  modes  comprising  surface 
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waves  in  both  the  lower  and  upper  half-spaces.  However,  for  <2  >  el>  M  COQ* 
sidered  here,  such  wave  modes  cannot  be  supported  and,  consistently  so,  Kj^n 
singularities  are  not  intercepted  in  the  deformation  of  the  SIP.  . 

The  expressions  pertaining  to  the  pole  contributions  discussed  above  are 


UU(T)  =  Re  { 


u  rk**  ..  -ik.TTfL  ^ — ,h[1}(Kuuup) 


JA  [ 

4*  Jk 


dkoe-^Hko)-  - - 


6k/^  *o) 


Af  =  integer  (15) 


H+mlu(t)  =  &e  S 


a  r 
4lr 


k,=KMVV  ) 
dk0  e~tk°r  I(k0) 


°°  •  KMLU  rxU) 


S  -  ki 


#1  (Kmlup) 


MLU 


*?=Kmlu  ) 


M  =  integer  (16) 


The  total  magnetic  field  may  then  be  evaluated  from  (10)  using  (8) -(9)  and 
(15) -(16).  This  is  what  we  shall  refer  to  as  the  tingle  deformation  approach. 
However,  the  double  integrals  in  (8)  and  (9)  are  not  easily  evaluated  due  to 
the  oscillatory  nature  of  the  integrand  as  a  function  of  real  k0.  The  double  de¬ 
formation  method  seeks  to  alleviate  this  problem  by  deforming  the  real  k0  path 
of  integration  to  the  positive  or  negative  imaginary  k0  axis.  At  the  same  time 
we  are  potentially  able  to  infer  physical  processes,  that  may  be  occurring,  from 
the  intermediate  results  obtainable  through  the  deformation.  We  reiterate  that 
the  deformation  in  the  k0  plane  is  performed  only  on  the  steepest  descent  path 
integrals  in  the  kp  plane  and  not  the  entire  original  field  expression.  As  such,  con¬ 
tributions  derived  from  various  singularities  in  the  k0  plane  should  be  interpreted 
as  contributing  to  the  time-response  of  these  integrals  which,  for  time-harmonic 
excitation,  correspond  to  both  the  direct  (from  saddle-point),  lateral  (from  branch 
cuts)  waves.  This  is  so  when  both  the  dipole  and  the  observation  point  are  on  the 
upper  boundary  in  which  case  the  saddle  point  coincides  with  the  branch  point 
at  kp  =  k0  [6]. 
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IV.  DEFORMATION  IN  THE  k0  PLANE 

Interchanging  the  order  of  integration  in  (8)  and  (9),  we  have 

■  S’-  j 

*!''«*•  “•’-'hliKSi’SKSill 


'*!-(* 7+W 


B*SDP2(t)  =  (  J  dq  J  dk0 e~»*«TJ(fe0)  ■  J~2  *  W) 

8ir  l  •'O  «/»  >/*2  -  W+  *«)2  ( 

•  j!"  (U,  + ) 

/<M(*2  +  *?»  *o)  ’•“(is  +  *9,  fco) J  J 

The  deformation  of  the  original  k0  integration  path  to  the  positive  or  negative 
imaginary  axes  is  determined  by  the  value  of  r  relative  to  p  and  also  to  Py/eii 
which  then  determines  the  convergence  of  the  resulting  deformed  integral. 


i 

t 

A 

l 

l 

»»  I 
~  1  i 

i 


•,  deformed  path 
r  <  P  ^ 


original  path  of 
°  integration  for  SDP1 


•  *,  deformed  path 

;  <•  <  ry/n  x 


i  |  original  path  of  » 

1 1  integration  for  SDP2  ' 


W 

v*i  *  1 


deformed  path 

'  r  >  „  y 


v/rj  -  l 


deformed  path 


Figure  4.  Original  and  deformed  paths  for  k0  integration  for  (a)  SDPl 
and  (b)  SDP2. 

We  find  that  [Fig.  4]: 

(i)  For  r  <  p,  the  k0  integrals  for  both  SDPl  and  SDP2  are  deformed  to  the 
positive  imaginary  axis  ( k0  =  iu,  u  >  0). 

(ii)  For  py/ti  >  t  >  p,  the  k0  integral  is  deformed  to  the  negative  imaginary  axis 
(ka  =  -iu,  u  >  0)  for  SDPl  and  to  the  positive  imaginary  axis  for  SDP2. 

(iii) For  r  >  the  k0  integrals  for  both  SDPl  and  SDP2  are  deformed  to 

the  negative  imaginary  axis  ( kc  —  -iu,  u  >  0). 
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Corresponding  to  (i).  ( ii )  and  (iii).  each  deformation  results  in  the  original,  real 
J bo  axis,  integration  path  being  replaced  by  the  positive  or  negative  imaginary 
axis  and  an  arc  of  infinite- radius  subtending  either  the  first  or  fourth  quadrant. 
Consistent  with  each  deformation,  the  integrations  over  the  arcs  result  in  zero  con¬ 
tribution  by  Jordan’s  lemma.  The  contributions  to  H0  SDP\lT)  an<*  ^<s> SDP^T) 
are  to  be  found  in  the  resulting  double  integral  over  q  and  u  and  the  residues 
of  any  singularities,  in  the  first  or  fourth  quadrants,  intercepted  in  the  deforma¬ 
tions.  These  singularities  are  characterized  by  the  zeros  of  g^ ( k0  *  iq,k0)  and 
4-  iq,k0)  for  SDP1  and  the  zeros  of  *<7i  fco)  and  rj^(k2  *9*ko) 

for  SDP2.  We  must  also  consider  the  singularities  associated  with  I(k0).  For 
SDPl,  there  are  branch  points,  corresponding  to  the  transformed  koz  and  kiz, 
at  k0  =  -iq/2,  and  at  k0  =  -iqjiy/t^  —  1)  and  iq/(^/i 2  -  1)  respectively.  The 
branch  cuts  chosen  are  shown  in  Fig.  4(a).  The  sign  of  the  square  roots  asso¬ 
ciated  with  the  transformed  koz  and  ib-2- ,  evaluated  anywhere  in  the  complex 
k0  plane,  should  be  determined  with  reference  to  the  fact  that  their  real  parts 
are  positive  on  the  original,  real-axis,  integration  path.  On  the  other  hand,  for 
SDP2,  the  branch  points  are  located  at  k0  =  ~iq/( 2y/(2  )>  corresponding  to  the 
transformed  k2z,  and  at  k0  =  —  * 9 / ( v/«2  f  1)  &nd  -l9/(\/*2  ~  1)  corresponding 
to  the  transformed  koz  ■  The  branch  cuts  are  shown  in  Fig.  4(b).  The  signs  of  the 
square  roots  associated  with  the  transformed  koz  and  k2t,  evaluated  anywhere 
in  the  complex  k0  plane,  are  again  determined  with  reference  to  the  fact  that 
their  real  parts  are  positive  on  the  real-axis. 

The  singularities  in  the  k0  plane,  determined  from  the  zeros  of  g^iko  +  iq,  k0) 
and  denoted  by  P\fg,  are  the  solutions  of 

In  Ro\R\2k  =k  *  ftdyjiT\kl  -  (k0  1-  i$)2  =  i(2M  +  1)7T  M  =  integer 

(19) 

The  zeros  of  /tM(fc0  —  ‘9-  k0),  denoted  by  P.vf/i  satisfy 

In  ^  -  i2dJtrlkl  -  {k0  +  iq)2  =  i(2M  +  l)ir  A/  =  integer  (20) 

.P01 

In  addition,  for  SDP2.  the  zeros  /tM(*2  +  iq,k0),  denoted  by  P.1//2  satisfy 

-i2 dJtrikl  -  (*2  +  iq)2  =  i(2.Vf  -K 1  )rr  .V/  =  integer  (21) 

and  the  zeros  r^"*(  denoted  by  P\fr  satisfy 

-  In  RaiRulk  =k2+w  *  *2  -  (k2  -  iq)2  =  i(2M  +  1  M  =  integer 

*  "  (22) 
where  all  square  roots  in  ( 19)  -  (22)  are  defined  to  have  positive  real  components. 
The  details  of  solving  (19) -(22)  are  examined  in  t9|. 

The  loci  of  P\fg  and  P\tf\  as  functions  of  qd  are  illustrated  in  Figs.  5(a)-(b) 
for  =  3.2  and  e2  =  80.  From  Fig.  5(a),  we  find  that,  for  each  M ,  the  locus 
for  P\jg  intersects,  if  at  all,  the  fc0  axis  at  two  frequencies  corresponding  to  k\f  \ 
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Figure  S.  Pole  lod  as  function  of  qd  in  k0d  plane.  No  pole  loci  are 
enclosed  for  Pjgfi  and  P\fr  in  the  deformation.  —  3.2,  e<i  = 

80.  (a)  Poles  P^g  for  M  =  0,  1,  2,  3.  (b)  Poles  P^f\  f°r 
M  =  1,  2,  3. 

and  ktfi-  The  locus  for  P\ff\,  on  the  other  hand,  intersects  the  k0  axis  at 
fcj/j.  We  denote  the  corresponding  q  values  by  qM\,  qm  and  qMS-  These 
critical  frequencies  are  identical  in  value  to  the  k^i,,kif2  and  kjy 3  in  (15)  and 
(16)  and  are  crucial  to  the  determination  of  causality. 

Neither  of  the  loci  for  Puf2  and  P&(r  exists  in  the  domain  of  deformation 
for  the  cases  we  consider  with  «2  >  <1  >  1-  As  such  they  do  not  contribute  for 
r  <  Ps/t2  and  do  not  pose  any  problems  in  terms  of  satisfying  causality. 
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We  find  in  particular  a  pole  locus  designated  by  M  =  0  for  P^ig  no 
correspondence  to  any  locus  intercepted  in  the  kp  plane  since  it  never  intersects 
the  real  k0  axis  and  lies  entirely  on  the  negative  imaginary  axis.  The  Appendix 
discusses  the  details  associated  with  the  existence  and  search  of  this  elusive  pole 
locus  and  its  contribution  to  the  magnetic  field.  Its  location  on  the  negative 
imaginary  k0  axis  implies  a  response  which  decays  exponentially  with  time  and 
hence  is  essentially  important  only  for  the  early-time  response.  The  corresponding 
results  and  implications  will  be  discussed  in  Section  VI. 

Causality  requires  that  the  total  field  solution  be  zero  before  the  time  it  takes 
the  electromagnetic  wave  to  travel  the  separation  distance  at  its  speed  in  the 
intervening  medium,  which  in  this  case  is  the  less  dense  of  the  two  interfacing 
media,  has  elapsed.  The  waves  excited  as  modes  and  attributed  to  the  kp  plane 
singularities  are  time-harmonic  modes,  each  of  which  is  valid  for  all  times  at  a 
single  frequency  and  are  therefore  non-causal.  The  summation  of  poles  on  the  kp 
plane  in  any  manner  also  does  not  result  in  a  causal  solution.  In  [4,9],  one  finds 
how  the  contributions  due  to  singularities  in  the  first  quadrant  of  the  k0  plane 
(intercepted  by  deformation  for  r  <  p)  will  combine  with  these  non-causal  modes 
to  give  a  zero  contribution  for  r  <  p.  The  fact  that  the  double  integrals  over  q 
and  u  for  the  corresponding  deformation  combine  to  give  a  zero  contribution  to 
the  magnetic  field  over  this  period  is  also  illustrated. 


IV.l  H^{r)  for  p^/ti  >  r  >  p 

With  causality  accounted  for,  we  next  consider  the  time  interval  Py/ei  >  r  > 
p.  Only  SDPl  is  deformed  downwards  while  SDP2  is  deformed  upwards.  The 
contributions  to  H^SDP\{T)  ^  com«  from:  B*Mg  B*Mfl  —  due  to 
P\fg  and  P\tf\\  H+iDP  —  from  the  singularities  due  to  the  source  function 
I{k0)i  HSDpyd  —  the  remaining  double  integral  over  q  and  u.  The  upward 
deformation  of  SDP2  results  in  Hsopy*  —  the  remaining  double  integral  over  q 
and  u,  i.e. 

B*SDPl(T)  =  ^2S*Mg  +  ^B+Mfl  +  B<plDP  +  H<t>SDP\*  (23) 

M  M 

H<t>SDPl(T)  =  B4>SDP7*  (24) 


where 


HpMg(T)  ~ 


-?*'< 

4w 


( P\fg  +  W)2 


yJPhg-^Hg  +  v)2 


rrm,,r  .  x  x  ft1  (Pm g  +  Pm g) 

■  B\  \{PMq  +  *g)p)-y--  l 


9k , 


(25) 
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]  'j  *  'Pm,X Tl{PMf\)~j  }PMfX  *  tq)— 

4*  I  ■'O  '/-Pi/1-(^M/1+*?)2 


•  B^\{PMn  +  iq)p)  g<M(P^^  +  Pyf- -  >  (26) 

£7-/«M(fco  +  *9»  ko) 

°k°  ko=Ptifl  . 

In  accordance  with  the  choice  of  source  function  (1),  /(ho)  has  a  pole  of  (n  +  1)  th 
order  in  the  fourth  quadrant  at  k0  =  kp,  then 

B*iDP(t)  =  |  iUm,  ^(k0  -  k,)"+'e-**rI(k 0) 

•  H  dq  -  Jrk%?--q)2-=H\l)((k0  +  »g)p) 

L  /“(ho  +  «g,  ho)  ftM(h0  +  ig,  h0)  1  ,27> 

[  StM(*o  +  >«,  h0)  /“(ho  +  tq,  h0)J  J 

The  expressions  for  the  remaining  double  integrals  take  the  form 

B4>SDP\*(t)  =  d(*  JQ 

■  ■..■i*-C^l-^s{l>(iU  -  «W  (28) 


The  expressions  for  the  remaining  double  integrals  take  the  form 


-  ur  -  « 


/“(fep,  -»«)  ftM(feP»  ~»u) 

9t*(kPi  ~*u)  ft*(kP'  “,u). 


' 

”  "£J**  ■/  *Jt 

k 


dxieUTI{iu) 


■  -!£±£ML  -g}1)(i(«  +  (29) 

y/W  +  «v^5)2  “  “2 

g«M(hp,  »*)  ^“(hp,  *«)  1 

.  /tM(fcP’*u)  rt*(hp»,«)Jjk#ssi(,+M>/5j)  J 

The  presence  of  the  Af  =  0  loci  for  Pj/,,  which  are  the  zeros  of  $<“(*(?  - 
u),  -iu),  on  the  negative  imaginary  axis  h0  =  — its,  modifies  the  corresponding 
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u  integration  in  (28)  to  be  a  principal  value  integral.  In  addition,  the  term  *2’  in 
square  brackets  in  (28)  may  be  dropped  in  the  computation  as  the  corresponding 
solution  is  purely  imaginary.  This  can  be  readily  shown  by  interchanging  the 
order  of  integration  and  performing  the  variable  transformation,  Q  =  q  -  u. 

To  summarize,  the  total  time-domain  solution  of  (5)  for  p^/ti  >  t  >  p  is 
expressible  in  terms  of  (15) -(16)  and  (23)-(29),  i.e. 

H+ip  <t  <  p^/ei)  =  ^  H*mvu(t)  +  51  S*MLU(T) 

M  id 

+  5Z  +  51  B*Mfl(T) 

id  M 

+  DP(T)  +  B4SDP\*(t)  +  B*SDP2*{t) 
where  all  components,  except  for  B^SDpid  and  SD P2* ,  comprise  single  in¬ 

tegrals  while  H+SDPld  and  HjSDPV  are  double  integrals  of  rapid  convergence 
by  construction. 


IV.2  B+(t)  for  r  >  p v/«£ 

We  now  consider  the  remaining  time  interval  r  >  py/l 3.  Both  SDP1  and 
SDP2  are  now  deformed  downwards.  The  contributions  to  ffysDPl(T)  will  come 
from:  H+  Mg  and  —  due  to  PMg  and  ~  from  the 

singularities  due  to  the  transformed  source  function  J(Jb<>);  and  BSDpxd  —  the 
remaining  double  integral  over  q  and  u.  The  expressions  for  these  contributions 
are  identical  to  those  in  (25) -(28). 

The  downward  deformation  of  SDP2  results  in  contributions  from  the  poles  of 
Pjd  f2  and  P\(r,  the  source  pole  at  k0  =  kp ,  and  the  contribution  from 
—  the  remaining  double  integral  over  q  and  u,  i.e. 

BtSDPliT)  =  5 Z,B<Hdg  +  5 +  B^XDP  +  B^SDPXd  (30) 
id  M 

Bj  SDPl(T)  =  5 Z,B*Mf7  +  5 2B<Hdr  +  B+IDP  +  H<f)SDP2  *  (31) 

id  id 


where 


I  i  M Pu/2)  /(PUPfr2+'^ 

"  I  «  y/Phf2-(PMflV£i+'l)2 


B\\Puf 2  +  'iM-0 


9tl(Pldf2y/yj  +  W»  PMf2) 


Qj^ft  (koy/Ui  +  W>  ko) 


ko-P\tfi 


(32) 
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r)  =  if  d,'-ir"'THl >«,)  (iWS?  +  <«>>  ■ 

4*  I  v/p^-(p«^+,5) 


rr(l)//«  ,  x  JtM(^Wrv^2  +  *9.  -PMr) 

•  ((•*  J/r V=*2  -«•  W)p)  g  - - - 

17  -M/ 


Qj~rt  (koy/Ui  +  *9.  *•) 


ko=Pi4r 


(33) 


We  have  stated  that  for  the  cases  under  study,  P\tf2  and  Pj^r  are  not  inter¬ 
cepted.  Equations  (32) -(33)  are  shown  here  for  the  sake  of  completeness  and 
are  applicable  in  cases  where  they  are  enclosed  in  the  deformation. 

For  the  source- pole  contribution,  we  have 


DP(.T) 


lim 


=  —  Re  {  i  lin 
4irn!  ( 

/, 


<T 

dJb» 


(kc  ~  kp)n+1e-i*»r/(Jfeo) 


dqM^akoy/s^  +  i^ph  (fc°^  +  w) 


+  if,  ka) 


\  ff{k0y/tr2^iq,k0) 

For  the  deformed  double  integral  for  SDP2, 


y/kj  -  (k0y/£i  +  W)2 

^(koy/i^  +  iq,  fceA  ) 
r?(koVUi  +  '<h  ko)J  J 


(34) 


y/W^VZi)2  ~  “2 


HJ  '(*(?  +  «v^2V) 


(35) 


To  summarize,  the  total  time-domain  solution  of  (5)  for  v  >  is  expressible 

in  terms  of  (15)-(16),  (23)  —  (28),  and  (32)-(35),  i.e. 

>  Pv/<2  )  =  £  +  Y  H*Mg(T) 

M  hi  M 
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+  Y  SdMf\(T)  +  BdlDp(T)  +  Y  S*Uf2(T)  +  Y  S+Mr(T) 

M  MM 

+  B*2  DPir)  +  B*SDpAt)  +  BdSDP2^T) 
where  again  all  components,  except  for  B^SDPld,  and  B^SDP2d,  consist  of 
single  integrals  and  B^sDP\*  BdSDP2 4  double  integrals  of  rapid  con¬ 

vergence  by  construction. 

V.  COMPUTATION  OP  POLE  CONTRIBUTIONS 

In  the  double  deformation  approach,  we  display  responses  corresponding  to  various 
sets  of  poles  denoted  by  the  integer  Af.  Each  Af  comprises  poles  from  the  kp 
and  k0  planes.  We  emphasize  that  the  k0  plane  poles  are  associated  with  the 
decomposition  of  the  double  integrals  characterizing  the  time  response  of  the  kp 
plane  steepest  descent  paths.  We  also  note  that  the  evaluation  of  the  source-pole 
contribution  does  not  fall  within  the  scheme  to  be  discussed  here  but  is  performed 
in  a  straightforward  manner  based  on  (27)  and  (34). 

For  r  >  p,  we  remark  on  the  scheme  of  evaluating  the  pole  contributions.  This 
is  the  only  range  of  time  to  be  considered  and  there  are  no  solutions  corresponding 
to  BM f 2  BMr  [Section  VI].  The  evaluation  of  the  pole  contributions  are 

greatly  eased  by  the  correspondence  between  positions  of  kp  poles  and  k0  poles 
at  certain  locations.  In  particular,  th*  real  frequencies  at  which  the  kp  plane 
pole  loci  intersect  the  steepest  descent  path  are  identical  to  the  intersections  of 
the  ka  plane  pole  loci  with  the  real  frequency  axis.  This  is  expected,  since  at 
these  points,  the  root  equations  (12)  and  (19)  are  identical  and  so  are  (13)  and 
(20). 

In  the  scheme  employed  [4,9],  for  each  Af,  we  locate  the  position  and  evaluate 
the  contribution  of  each  pole  of  the  locus  on  the  kp  as  well  as  the  kQ  planes. 
Refering  to  Fig.  6,  which  depicts  the  situation  for  a  general  Af,  we  begin  the 
search  for  poles  in  the  k0  plane  starting  with  the  Pjjg  branch  for  qd  =  0.  At 
q  =  9Af2i  the  ko  pole  will  be  fcj/j  corresponding  precisely  to  the  value  for  k0 
when  the  kp  plane  poles  for  the  Af  th  locus  ‘enters’  the  ‘region  of  interest’  [Fig. 
3].  Moving  then  to  the  kp  plane  and  continuing  on  the  corresponding  kp  pole 
(K\fuu)  locus  [(15)]  by  lowering  the  frequency,  we  eventually  leave  the  region 
of  interest  at  Kmju  =  km  +iqMl-  *Afl  corresponds  to  the  frequency  at  which 
the  Pf4g  pole  lod  ‘returns’  to  the  fourth  k0  quadrant  (when  q  =  )  for 

consideration  when  r  >  p. 

The  computation  continues  on  the  locus  [(19)]  down  into  the  fourth 

quadrant  and  converges  quickly  since  the  integrand  has  a  e~4^  dependence  and 
is  also  proportional  to  ck»(r~p)  for  negative  kp. 

We  return  to  the  point  for  qd  —  0  for  the  Pm/1  branch  (which,  by  proper 
choice  of  Af ,  may  be  coincident  with  the  starting  point  for  the  P^g  branch)  and 
repeat  the  same  procedure  with  the  Pj^ ^  [  (20)  ]  branch  and  locus  [  (16)  ] 

instead.  For  the  Af  th  locus,  therefore,  the  computation  terminates  with  the  in¬ 
tegration  over  the  K^LU  P°le  contributions  in  the  k0  plane.  These  poles  lie 
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Figure  6.  Example  of  summation  path  of  pole  contributions  corresponding 
to  general  Af. 

closer  to  the  real  kp  axis  and  convergence  of  the  pole  contributions  is  slower  than 
for  those  in  the  ka  plane. 

Figure  6  shows  the  general  situation.  From  the  nature  of  the  proof  of  causality 
[4,9],  we  find  each  Af  response  to  be  causal  if  computed  in  accordance  with  the 
figure.  In  certain  cases,  for  low  orders  of  Af,  the  Pj^g  branch  may  not  intersect 
the  inversion  contour  on  the  k0  plane,  implying  that  no  corresponding 
poles  are  in  the  region  of  interest  and  thus  necessitating  a  slight  modification  to 
(19).  This  feature  can  easily  be  taken  into  account.  In  particular,  the  special  case 
of  Af  =  0  is  considered  in  the  Appendix. 


VI.  RESULTS  AND  DISCUSSION 

In  Figs.  7-11,  we  present  results  obtained  for  the  time-domain  magnetic  field  from 
the  double  deformation  method  in  comparison  with  existing  methods  as  well  as  for 
various  sets  of  parameters.  Unless  explicitly  stated  otherwise,  we  shall  consider 
the  source  function  of  the  form, 

*(<)  “  /o<sin(wo0«~ao<u-l(0  (36) 

The  initial  response  due  to  such  a  function  should  begin  from  a  zero  value  [1],  so 
that  a  convenient  check  is  immediately  available  for  our  numerical  results.  I{k0) 
has  a  pair  of  second  order  poles  at  k0  =  ±1.0  —  »0.5,  of  which  only  that  in  the 
fourth  quadrant  is  of  interest  and  will  constitute  the  source-pole.  To  allow  ready 
comparisons  with  alternate  methods,  we  consider  the  lossless  case.  In  addition, 
the  parameters  used  shall  be  =  3.2,  <2  =  80,  u»o  =  1.0c,  ao  =  0.5c,  d  =  1  m 
and  p  =  3  m  for  most  of  the  cases  illustrated. 

In  Fig.  7,  we  show  the  results  obtained  via  the  single  deformation  method. 
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Figure  7.  Time-domain  magnetic  field  versus  r  =  ct  for  VED  on  two- 
layer  medium  using  single  deformation  approach.  t\  =  3.2, 

«2  =  80,  p  =  3  m,  d  =  1  m,  i(t)  =  70tsin(ct)  • 

(*)  poles  contribution,  (b)  SDPl  and 

SDP2  contribution. 

Figures  7(a)  and  (b)  respectively  depict  the  response  due  to  the  kp  poles  and  the 
steepest  descent  path  integrals.  The  two  responses  are  comparable  in  magnitude. 
The  sum  or  total  time  response  is  plotted  in  Fig.  9(a)  in  comparison  with  the 
double  deformation  result. 

Figurer  8(a)-(f)  show  the  various  contributions  to  the  double  deformation  re¬ 
sults.  Figure  8(a)  shows  the  source-pole  contribution  which  follows  the  oscillating 
and  decaying  nature  of  the  source  function.  Figures  8(b)-(f)  show  the  responses 
corresponding  to  the  various  pole  loci  M  and  evaluated  according  to  the  scheme 
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Figure  8.  Magnetic  field  contributions  versus  r  =  ct  for  VED  on  two-layer 
medium  using  double  deformation  approach.  el  =  3.2,  e2  = 

80,  p  =  3  m,  d  =  1  m,  t(f)  =  fo<sin(ct)e-0'5c<u_i(t).  (a) 
Source-pole  contribution,  (b)  Lowest  order,  Af  =  0,  poles  con¬ 
tribution. 

described  in  Section  V.  In  particular,  we  note  that  the  Af  =  0  set  in  Fig.  8(b)  is 
essentially  an  early  time  response  due  to  its  location  on  the  negative  imaginary 
k0  axis.  We  also  find  that  the  contributions  fall  rapidly  with  increasing  Af,  and 
Af  <  5  should  be  sufficient  to  achieve  reasonable  convergence  in  results.  The 
total  response  after  summation  is  shown  in  Fig.  9(a). 

We  have  not  shown  the  double  integral  results  corresponding  to  sdpi4 
H<isDP2u  d  because  they  are  found  to  be  negligibly  small  for  all  cases  investigated 
and  appear  unnecessary  for  reasonable  accuracy.  Moreover,  in  [9],  it  is  argued  that 
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Figure  8.  (cont’d) — Magnetic  field  contributions  versus  r  =  ct  for  VED 

on  two-layer  medium  using  double  deformation  approach,  t\  = 
3.2,  <2  =  80,  p  =  3  m,  d  =  1  m,  i(t)  =  /otsin(c<)e-0'5ctu_i(t). 

(c)  AT  =  1,  (d)  M  =  2,  (e)  M  =  3,  and  (f)  M  =  4. 

the  double  integrals  evaluate  to  zero  for  the  non-dispersive  half-space  problem. 
This  implies  that  ov'  neglect  of  the  small  contribution  of  the  double  integrals 
becomes  increasingly  justified  for  larger  layer  thickness  d  in  the  lossless  case. 

The  accuracy  cf  the  double  deformation  results  is  verified  by  comparison  with 
the  single  deformation  waveform  in  Fig.  9(a)  and  with  the  explicit  inversion 
method  [1]  in  Fig.  9(b).  Figure  9(a)  confirms  the  validity  of  the  double  defor¬ 
mation  results.  Further  confirmation  is  obtained  by  matching  with  the  half-space 
results  [9]  over  an  initial  time-window  and  to  be  shown  in  Fig.  10(a).  The  com¬ 
parison  with  the  explicit  inversion  results  shows  slight  discrepancies  in  value  if 
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Figure  9.  Total  time-domain  magnetic  field  versus  r  =  ct  for  VED  on 
two-layer  medium.  «i  =  3.2,  <3  =  80,  p  —  3  m,  d  =  1  m, 
*(<)  =  Jot»in(ct)e-0,5c<u_i(t).  (a)  Single  deformation  and  dou¬ 
ble  deformation,  (b)  Explicit  inversion  and  double  deformation. 

not  in  trend.  In  the  explicit  inversion  scheme,  the  step  response  is  first  computed 
and  subsequently  convolved  with  the  source  function  to  obtain  the  required  re¬ 
sponse.  Much  relies  on  the  accuracy  of  the  step  response  which  has  to  be  evaluated 
from  an  integral  expression  [1].  In  Fig.  10  we  compare  the  half-space  solution, 
which  may  be  readily  obtained  from  the  closed-form  impulse  response  [10],  with 
the  double  deformation  results  for  d  =  1  m  and  d  =  0.2  m.  The  initial  deviation 
from  the  half  space  solution  indicates  the  first  arrival  due  to  reflection  off  the 
lower  boundary. 


154 


Poh  and  Kong 


Figure  10.  Magnetic  field  contributions  versus  r  =  ct  for  VED  on  two- 
layer  medium  (double  deformation)  and  a  half-space,  d  — ♦ 
oo  (convolution  with  closed-form  impulse  response).  = 
3.2,  =  80,  p  =  3  m,  i(t)  =  /otsin(cf)e-0-5c<u_i(f).  (a) 

d  =  1  m,  (b)  d  =  0.2  m. 

In  Fig.  11(a)  we  show  the  results  for  single  and  double  deformation  for  the 
case  of  <2  =  12,  all  other  parameters  unchanged.  Due  to  the  smaller  dielectric 
contrast  at  the  lower  boundary,  the  effect  of  the  first  reflected  wave  is  not  as 
pronounced  as  in  the  previous  case  of  «2  =  80.  In  addition,  we  have  again  not 
included  any  of  the  double  integrals  corresponding  to  §DP\ 4  SD P2u  d 

in  our  computation.  We  find  that,  in  comparison  with  single  deformation  results, 
this  does  not  affect  the  accuracy  of  the  double  deformation  method. 


Sr-f 

f. 
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Figure  11.  Magnetic  field  contributioiu  versus  r  =  ct  for  VED  on  two- 
layer  medium,  p  =  3  m,  d  =  1  m,  t\  =  3.2.  (a)  «2  = 

12,  t(f)  =  /0 tsin(ci)e-0-5c,u_i(f).  Single  and  double  defor¬ 
mation.  (b)  «2  =  80,  *(<)  =  Io  tn  sin(c<)e-0-5rtu_i(<),  with 
n  =  1  and  2.  Double  deformation. 

The  effect  of  using  a  smoother  source  function  i(t)  =  J0<2  sin(c<)e-0'5elu_i(f) 
is  illustrated  in  Fig.  11(b),  in  comparison  with  the  source  described  by  (36).  Here 
the  most  outstanding  difference  is  the  slower  and  smaller  change  in  the  waveform 
due  to  the  first  reflected  arrival.  The  smoothness  of  the  source  accounts  for  the 
delayed  effect  and  also  allows  the  falling  portion  of  the  half-space  solution  to 
smooth  out  much  of  the  early  contribution  of  the  reflected  signal. 

We  have  illustrated,  by  means  of  a  few  examples,  the  feasibility  of  the  dou¬ 
ble  deformation  method  applied  to  transient  dipole  radiation  over  a  two-layer 
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medium.  In  comparison  with  a  previous  attempt  in  this  direction  [4],  the  identi¬ 
fication  of  the  M  =  0  pole  locus  enables  us  to  obtain  a  complete  response  to  the 
time-domain  problem. 

APPENDIX:  EXISTENCE  OP  P0f  LOCUS 

In  this  appendix,  we  show  how  the  pole  locus,  designated  M  =  0  for  PMg, 

satisfying  the  equation  jj*  =  0,  on  SDPl,  is  found  to  exist  on  the  negative 
imaginary  k0  axis. 

For  the  sake  of  illustration,  we  consider  (2  >  ei-  la  the  k0  plane,  for  SDPl, 
koz,  k\x  and  fcjx  may  be  written  as 


so  that  there  are  branch  points  at  k0  =  —iq/2  corresponding  to  kox  and  ka  = 
i?/(\/«2  ~  i)  and  ko  =  -iq/{y/ii  +1)  corresponding  to  *2*-  The  locations  of 
these  singularities,  the  corresponding  choice  of  branch  cuts  and  the  critical  points 
at  ka  =  iq/{y/e\  -  1)  and  k0  =  -iq/( y/t[  +  1)  for  k\t  are  shown  in  Fig.  A.l. 

Refering  to  Fig.  A.l,  we  may  divide  the  negative  imaginary  axis  into  four 
intervals  namely,  0  <  X  <  q/(v^+l),  q/{y/i 2  +  1)  <  X  <  q/{y/ti+l),  g/(v/«T+ 
1)  <  X  <  qj 2  and  X  >  qj 2.  We  have  let  k0  =  —iX,  where  X  >  0.  We  note 
that  only  in  the  interval  X  >  q/2,  ait  all  of  kot,  k\z  and  purely  imaginary, 
i.e. 


so  that  all  of  J2oi*  ^12  e*2*1^  real. 

We  find  that  over  this  portion  of  the  imaginary  axis,  gf*  =  0  and  fj*  =  0 
become  real  equations  in  X,  i.e. 

gt^iko  +  iq,k0)  =  \  +  Ro\Ri2ei2kltd  =0  (A.l) 

k0=ko+iq 

/,M(*o  +  ■»,*«)  =  *,,+«!,  =0  (A. 2) 

kp=k0+tq 

are  real  equations,  for  k0  =  -iX  and  X  >  q/2.  Consequently,  (A.l)  and  (A. 2) 
may  be  readily  solved  using  algorithms  for  real  functions.  This  would  be  far  sim¬ 
pler  than  attempting  to  determine  the  pole  positions  from  the  complex  equations 
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as  given  in  (19)  and  (20).  The  real  solutions,  multiplied  by  — would  give  P0? 
and  PQf  (M  =  0)  respectively.  We  find  that  only  (A.l)  supports  any  solution 
and  that  there  are  two  sets  of  solution  for  varying  q,  one  lying  entirely  between 
X  =  q/ 2  and  X  —  q  and  the  other  entirely  in  the  range  X  >  q.  Refering  to  (25) 
we  find  that  only  the  second  set  gives  a  non-zero  contribution  to  the  real  magnetic 
field  and  it  is  this  set  that  we  refer  to  as  Pq9  and  is  depicted  as  the  M  =  0  loci 
in  Fig.  5(a). 
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Figure  A.l.  Branch  points  and  branch  cuts  in  the  kQ  plane  for  SDP1. 


The  magnetic  field,  as  a  modification  to  (25),  is  therefore 


iWr  >  *>  *  -  £*•  |  if0  dqe~iPM'rr{PM')ffif: 


{Pm g  +  w)2 


(PMo  +  *9)2 


i  ftA{PM9+i<I>  PM9) 


dk. 


■9t(ko  +  iq ,  *o) 


ko-PMg  )  M-q 


where  analytical  simplification  may  be  performed  before  computation  by  recog¬ 
nizing  that  P0g  is  purely  imaginary  and  that  the  loci  lies  entirely  on  k0  —  —iX, 
and  X  >  q.  Half  the  residue  is  taken  due  to  the  location  of  the  poles  and  the 
nature  of  the  deformation. 
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Time  Domain  Analysis  of  Nonuniformly  Coupled  Line  Systems 
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Massachusetts  Institute  of  Technology 
Cambridge,  MA  02139,  USA 

Abstract-  A  general  method  of  analyzing  the  time- domain  bi-directional  coupling  of  a  pair 
of  nonuniformly  coupled  dispersionless  transmission  lines  is  presented.  The  transmission 
line  equations  are  decoupled  using  the  method  of  characteristics  and  the  equations  are 
solved  iteratively.  In  the  cases  with  linear  loads,  the  unit-step  response  can  be  obtained 
in  closed-form  to  the  first  order,  and  arbitrary  excitations  can  be  handled  by  convolution. 
Numerical  integration  for  the  cases  with  nonlinear  loads  is  also  shown  to  be  more  efficient 
than  integration  based  on  the  original  coupled  partial  differential  equations. 

INTRODUCTION 

Previous  analyses  on  coupled  line  systems  have  often  emphasized  mono-directional 
(codirectional  or  contradirectional)  coupling  problems  [l]-[5].  The  interest  in  such 
problems  is  due  to  the  fact  that  the  mono-directionally  coupled  line  system  is  a 
very  accurate  model  for  many  practical  situations  such  as  coupled  TEM  trans¬ 
mission  lines  and  dielectric  waveguides,  and  that  the  model  simplifies  the  analysis 
of  the  corresponding  coupling  problem.  With  increasingly  compact  packaging  of 
electronic  circuits,  especially  in  the  high-speed  digital  integrated  circuits  used  in 
computers,  crosstalk  between  adjacent  signal  lines  in  circuits  becomes  a  significant 
problem  in  that  it  limits  the  speed  of  circuits  as  well  as  the  maximum  packaging 
density.  It  is  well  known  that  crosstalk  is  caused  by  electromagnetic  coupling. 
However,  the  coupled  system  consisting  of  signal  lines,  i.e.  interconnecting  lines 
in  integrated  circuits,  is  no  longer  a  mono-directionally  coupled  system  as  a  result 
of  the  inhomogeneities  caused  by  the  complex  configurations  [6]- [8]. 

The  crosstalk  in  low  speed  integrated  circuits  or  low  frequency  printed  circuit 
boards  is  usually  dealt  with  using  a  lumped  parameter  model  [9j-[10].  Yet  recent 
progresses  in  digital  integrated  circuits  have  extended  clock  speeds  into  the  3  GHz 
region,  with  rise  times  of  signals  dropping  below  200  ps,  hence  the  bandwidths 
of  the  signals  on  these  circuits  have  increased  to  several  GHz  or  more.  Obviously 
lumped  parameter  analysis  is  no  longer  accurate  and  all  interconnecting  lines  in 
circuits  must  be  treated  as  microwave  transmission  lines.  To  be  specific,  they 
should  be  modeled  as  nonuniformly  coupled  transmission  line  systems,  and  the 
coupling  is  bi-directional.  This  is  harder  to  deal  with  than  the  mono- directional 
coupling  problem.  In  [12j-[14]  the  authors  used  approximate  differential  equa¬ 
tions  under  the  assumption  of  weak  coupling  to  discuss  transmission  properties 
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and  coupling  effects  of  high-speed  integrated  circuits  and  printed  circuit  boards. 
With  the  method  of  characteristics,  computer  simulation  based  on  the  equivalent 
circuit  of  individual  modes  is  commonly  employed  to  analyze  transients  in  coupled 
transmission  line  systems,  and  gives  rather  accurate  results  [15j-[18j,  but  such  a 
technique  could  be  time-consuming  when  nonuniform  coupling  is  involved.  In  fact, 
until  recently  most  analyses  of  the  transient  responses  in  integrated  circuits  were 
based  on  uniformly  coupled  transmission  line  models.  Yang  et  al.  [19]  solved  the 
nonumformly  coupled  transmission  lines  problem  with  an  elegant  time  domain  for¬ 
mulation,  but  only  the  contradirectional  coupling  can  be  handled  without  much 
effort.  In  addition,  some  special  kinds  of  transmission  line  systems  with  peri¬ 
odically  varying  coupling  properties  were  analyzed  with  the  odd  and  even-mode 
equivalent  circuits  and  transform  technique  [20]- [21]. 

The  purpose  of  this  paper  is  to  generalize  the  transient  analysis  method  given  in 
[19].  A  new  set  of  variables  is  introduced  such  that  the  transformed  equations  are 
simpler  and  both  the  codirectional  and  contradirectional  coupling  can  be  easily 
calculated.  By  virtue  of  the  formulation,  causality  is  preserved  and  each  higher 
order  term  in  the  perturbations!  series  can  be  interpreted  as  a  partial  reflection 
along  the  lines  due  to  nonuniform  coupling.  Since,  in  practice,  the  terminations 
for  interconnecting  lines  in  integrated  circuits  can  be  linear  or  nonlinear,  we  shall 
investigate  both  cases.  Solutions  accurate  to  the  first  order  in  spatial  derivatives 
of  the  coupling  coefficients,  which  is  analytically  manageable,  will  be  presented. 


FORMULATION  AND  CHARACTERIZATION 


Consider  two  nonuniformly  coupled  non-dispersive  transmission  lines  as  depicted 
in  Fig.  1,  where  af(z,t)  and  a^{z,  t)  are  the  forward  and  backward  modes  on 
each  of  lines  1  and  2  [22].  Define  the  column  matrix  A  as 

A  =  (aJ-,af,o^,o^)‘ 


then  the  general  coupled  mode  equation  can  be  written  in  matrix  form 


1 

0 

~k+ 

dA 
dT  + 

'  0 

Pi 
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0  ' 
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-*+ 

-1 

-k- 

Jfe- 

1 

*+ 
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Pi 

0 

0 

0 

0 

0 

0 

P2 

*_ 

■ 

0 

-1 

0 

0 

P2 

0 

dA  1 
dz  v 


where  we  assume  that  two  lines  have  identical  phase  velocity 

1  1 


v  = 


>1  =  0 


(1) 

(2) 


and  that  the  propagating  modes  on  each  line  are  quasi-TEM.  L{  and  C{  (*  =  1,2) 
are  the  inductance  and  capacitance  per  unit  length  of  each  line,  and  are  all  func¬ 
tions  of  position  z .  The  elements  fc+  and  fc_  in  the  coupling  coefficient  matrix 
are  referred  to  as  the  codirectional  and  contradirectional  coupling  coefficients  of 
the  coupled  transmission  lines  respectively.  They  are  related  to  the  capacitive 
coupling  coefficient  kc  and  inductive  coupling  coefficient  ki  via 
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Figure  1.  Two  nonuniformly  coupled  transmission  lines. 


k+  = 


*c(*)  -  M*) 


k  _  M*)  +  kL(z) 
2 


(3) 


and  it  is  well  known  that  k c  and  ki  are  associated  with  the  mutual  capacitance 
C12  and  the  mutual  inductance  £12  in  the  following  forms 

Cn(z) 


kc  = 


y/CMCfc)  ’ 


kL  =  -M£L 


w 


s/Ll MiiW 

The  elements  pi  (i  =  1,  2)  of  the  coefficient  matrix  account  for  the  effects  of 
nonuniform  coupling: 

1/2 


—  k+  —  k- 
TkZ 


(*  =  1,2)  (5) 


where  Zx(z )  =  Lx(z)lCx(z)  and  Zqx  are  the  self-impedances  of  the  same  trans¬ 
mission  lines  with  infinite  separation,  i.e.,  lines  without  coupling. 

The  matrices  in  (1)  can  be  recast  in  partitioned  form 


s+ 


* 

0* 

_ -1 

9A 

>1 

:  0 

•  •  •  •••  ••• 

dt  + 

•  •  • 

•  •  •  •  •  • 

L  Ko  :  K ,. 

.  0 

:  Pi- 

A  =  0 


(6) 


where  the  coupling  coefficient  matrix  K  consists  of  the  following  elementary 
blocks: 


Kv  =  - 
w 


1  0 
0  -1 


and 


and  the  submatrices  Pi  and  P2  are 

0 


*1  = 


l  Pi 


Pi 

0 


and 


*o  =  - 
V 


P2  = 


—k+  -k- 
k-  k+ 


0 

P2 


P2 

0 


(7) 


(8) 
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The  coupled  mode  equation  (1)  can  be  cast  into  normal  form  by  diagonal¬ 
izing  the  coupling  coefficient  matrix  K  in  two  steps.  First,  using  the  rotation 
transformation 


2o  ^  ...  ... 


where  I  is  the  identity  matrix,  we  eliminate  the  off-diagonal  submatrices  in  K 


T0AT0-1  = 


K9  +  K0 


L  0  :  K,-K0  J  LO  :  K2  J 

Next,  we  diagonalize  the  submatrices  K\  and  K2  in  (10).  This  can  be  performed 
by  using  the  following  transformation  matrices 
1 

2^(1 -m2-h 

’ V(i -t- V(i-t+)rir  -*+)  +  *-' 

•y/(l  —  k+)  —  k-  —  ^/(l  —  k+)  -(-  k-  y/(l  —  k+)  +  fc-  +  \/(i  ~  ^+)^“  k- 


2{/(i  +  *+)2-*l 

\/(l  4-  fc+)  ■+•  fc—  +  \/(l  +  fc+)  —  k—  \/(l  +  k+)  +  kZ  —  y/(1  +  fc+)  —  k— 
\/(l  +  fc+)  +  k-  —  y/(l  +  &+)  ~  k—  \/(l  +  k+)  -I-  k-  +  y^(l  +  k+)  —  k- 


The  corresponding  diagonalized  submatrices  will  be 


\  ®  =  XKT~ 1 

0  XT  1  li 


_  I  [  [(l=F*+)3-*-] 


o  - [,i t*+>2 - *i]i/2J ( 

(*  =  1,2) 

Here  the  top  sign  is  for  »  =  1  and  the  bottom  sign  for  *  =  2 .  The  composite 
transformation  matrix  T  which  reduces  the  coupling  coefficient  matrix  K  into 
diagonal  form  is  then 


(14) 
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The  diagonal  matrix  A  is  related  to  K  by 

Ai  i  0 
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=  TKT~l 
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(15) 


0  :  A2  J 

Applying  transformation  T  to  (1),  and  after  some  algebraic  manipulations,  (1) 


becomes 


dV 
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0  ;  a2J 
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dv 

dt 


+  P2)T~l 
T2(Pl-P2)Tl- 1 


-faiPi-FtWi 


-1 
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-1 


V  (16) 


where 

(17) 

When  both  of  the  coupled  transmission  lines  have  the  same  parameters,  (16)  is 
simplified  to 
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*  *2  (l:pfc+)2-*2  2  (* -1’2) 


(18) 

(19) 

(20) 


Equation  (18)  is  effectively  decoupled.  The  energy  exchange  between  the  forward 
wave  V+  and  backward  wave  occurs  only  because  of  nonuniform  coupling.  If 
the  impedance  Zqi  of  line  1  is  independent  of  position  z  ,  b\  in  (18)  will  be  zero, 
and  the  forward  wave  Vf  and  the  backward  wave  V^~  are  completely  decoupled. 

Equation  (18)  can  be  further  simplified  by  introducing  four  families  of  charac¬ 
teristic  curves  which  satisfy  the  following  ordinary  differential  equations  [23] 


and 


dz 


(*  =  1,2) 

The  explicit  expressions  of  the  characteristic  curves  are 


—  =  \- 
dz  * 


(21) 


C,+  :  t  —  J  dz  =  constant 


and 
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-/**■- 


constant  (*  =  1,2) 


(22) 


On  these  four  characteristic  curve  families,  the  corresponding  equations  of  (18) 
become 

dV+ 

17  *  iiV‘~ 

and 


dVr 

dz 


=  6iV.+  (i  =  1,  2) 


(23) 


It  is  known  from  the  expressions  of  the  characteristic  curves  that  the  solutions  of 
V+  and  V~  represent  wave  travelling  in  the  positive  and  negative  z  -directions 
respectively.  Equation  (23)  further  demonstrates  the  fact  that  there  is  energy 
exchange  between  positive-going  and  negative-going  waves  between  two  bound¬ 
aries,  or  we  ran  interpret  it  as  continuous  reflection  along  the  lines.  This  type  of 
reflection  is  due  to  nonuniform  coupling. 

The  solution  to  (23)  are  to  be  constructed  iteratively,  making  use  of  time¬ 
marching  properties  of  wave  equations  [23].  Assuming  that  the  solutions  before 
time  t  have  been  obtained  at  any  position  along  the  lines,  then 

(24«) 

Vj±W(M)  =  >'S(*.0.!i<l>+  /W>  (246) 

Jzi0  Jz' 

(i  =  l,2;  n  =  0, 1, 2, 3, . . .) 

where  stand  for  solutions  to  (23)  with  the  right-hand  sides  equal  to  zero, 
which  are  simply  uniform  travelling  waves,  and  z^q  and  =  t  —  J^Q  dz'  A*  can 

be  chosen  such  that  f*o  <  *  •  Both  V-q (ziQ>*«o)  and  the  integrand  in  (246) ,  which 
is  taken  along  the  characteristic  curve  between  (zio,fio)  aad  (z,t),  are  obtained 
before  time  t .  Therefore  calculating  the  solutions  to  (23)  becomes  very  straight¬ 
forward.  Note  that  the  notations  stand  for  the  n  th  order  approximation, 

rather  than  the  nth  order  term  in  the  perturbational  series,  as  used  in  [19].  Since 
6j  and  6'2  are  proportional  to  Jfe+  and  k!_  ,  the  iteration  scheme  converges  to 
the  correct  answer  provided  that  kc  and  k i  are  slowly  varying  continuous  func¬ 
tions  of  Z .  The  general  discussion  on  convergence  can  be  found  in  [23],  Chapter 
II,  §  3.  The  solutions  for  forward  modes  at  and  backward  modes  a”  on  coupled 
transmission  lines  are  obtained  by  applying  the  inverse  transformation  of  (17) 

A  =  (a+,af,a2+,a2-)<=T-1K  (25) 


SOLUTION  FOR  THE  CASE  OF  LINEAR  TERMINAL  CONDITIONS 

We  shall  discuss  the  solution  of  two  transmission  lines  possessing  identical  pa¬ 
rameters  and  terminated  with  matched  loads  at  their  ends,  i.e.  Zt  =  Z( 0)  and 
Zi  =  Z(t),  and  excited  by  a  step  voltage  source  /(f)  =  u(t) .  According  to  the 
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relations  between  af 

and  voltage  ej  and  current  *,-  on  the  lines  [22] 

(26) 

The  initial  conditions 

are  zero,  and  the  boundary  conditions  will  be 

a+(0,f)  =  u(<) 

(27a) 

af(f,<)  =  0 

(276) 

a+(0,t)  =  0 

(27c) 

and 

«£W)  =  0 

(27d) 

The  boundary  conditions  of  variables  v*  can  be  obtained  from  the  transforma¬ 
tion  relation  as  expressed  in  (17).  For  convenience  of  expressing  the  boundary 
conditions,  we  rewrite  (17)  in  the  following  form 


(28a) 

(286) 
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From  (27)  and  (28),  the  boundary  conditions  of  V  are 
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V2-(f,t)  =  f22(f)V+(/,0 


where 


*i(*)« 


(<ia(«))< 


===tt(f)  +  il2(0)V2o(0,t)  (30c) 
+  *- 

(30d) 

(i  =  1,  2)  (31) 


110 


Gu  et  al. 


are  the  reflection  coefficients  at  the  boundaries. 

Based  on  (24),  we  can  use  a  similar  method  given  in  [19]  to  obtain  closed 
form  solutions  to  V  and  A  under  conditions  (30).  The  characteristic  curves  Cf 
associated  with  V*  are  straight  lines  as  in  [19].  Since 

A*  =  ±—  (32) 

«0 

where  v<)  is  the  phase  velocity  of  the  waves  propagating  along  the  uncoupled 
transmission  line,  the  corresponding  characteristic  curves  are 

t  ^  t\{z)  =  constant  (33) 


where 

Tl(z)  =  - 

But  the  characteristic  curves  C *  are  usually  not  straight  lines,  they  are 

t  ^  T2(*)  =  constant  (34) 


where 

*  r i  (i + k+(z))2 - ki(z') 

Z)'h  d  «0  ^  (1  -  k+(z'))2  -  fc2(z') 


(35) 


Along  these  characteristic  curves,  the  zeroth  order  solutions  to  V ,  which  corre¬ 
spond  to  the  solutions  of  differential  equations  (23)  with  zero  right-hand  sides, 
jure  waves  of  constant  junplitudes, 

„+(»),.  _ 2^(1 -M»»8- *-(»)» 

*•”  '  ’  ’  MO) +  *-(«)  +  -  Mo)  -  *-(0) 

(<  -  i  -  2>r,) 

„  _  »Si(<)t'(»-M0)),-*-(Q? 

1  '  yi  - 1+( o)  +  mo)  +  \fi  -  fc+<0)  -  mo) 

£>1(0)fi1(<)|-’u  (l  +  i  -  2(j  +  1)T,) 

„+(»),. _  2y(i  +  Mo));-t-(o); 

2'*  1  ’  '  v/1  +  MO)  -  MO)  +  Vl  +  *+(0)  +  *_(0) 

n 

£[J*2(°)/?2(/)]M*  -  r2(z)  -  2jTt) 

j= 0 

and 

„  _  2j?2(oy(i  +  M0))2-fc-(Q)it 

2’n  l2’  ’  y/l  +  *+(0)  -  *_(0)  +  y/l  +  Jb+(0)  +  M 0) 


£(*2(0)*2(*)]>«(t  -  r2(z)  -  2 (j  +  1)T2) 
j=o 


i 


* 
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(36) 

where 

Tl  =  f 

v0  _ 

r  =  tldz' -  /(1  +  M«'))2 -*!(**) 

2  Jo  «0  V  (!  -  M*'))2  -  *2  (*') 

and  t  is  the  length  of  the  coupled  region. 

Higher  order  solutions  to  V  can  be  obtained  by  making  use  of  the  iterative 
formula  (24).  The  general  form  of  the  first  order  solutions,  which  are  directly 
derivable  from  the  zeroth  order  solutions  y^0) ,  are  as  follows: 


In  region  (  nf  ):  2 nl*  -f  r* (z)  <  t  <  2(n  4-  1)2*  -  t*(z) 

%%  0  -  K.(0> *>,.<(*■  <.)  +  Vyi '2, 

n ;(«)) 

In  region  (  nb)  :  2(n  4-  1)T*  -  r*(z)  <  t  <  2(n  4-  1)T*  4-  r*(z) 

»S*W) = + vi;n,/1,K„,«-n) 

+ '•(')  - r') 

(<-!.*)  (37) 

and  y  ±(1)  =  0  in  the  region  below  the  characteristics  passing  through  (z“0,t“0) . 
Here 

if  1  -F  ^+(z)  T  k-(z)  1  T  k+(zi)  i  k-(zi) 

Rk*,iiz  'z*)=4|lni 


+ 

4-  In 


■f  ^>+(z)  i  (z)  1  T  ^+(z»)  T  k— (z») 

I  1  ~  *+(«)  +  M«) 1  ~  Mzi)  ~  M*»)  1 

"l  -  fc+(z)  -  fc+(z)  1  -  M*«)  +  M*»)  J 

Z0(z) 


Z<>(zi) 


(i  =  1,2) 


(38) 


z+  =  r'1 

z«,n  r» 


=  rr1 


jWz)  +  («-2;r,)) 


»  =  o(*  +  r<(z)]  +  nT, 


(39a) 


t{n(z)-{t-  20  +  1)7*)) 


,  t-|l  =  5[«-'ri(z)]  +  (n  +  1)r*  (396) 


where  r”1^)  is  the  inverse  function  of  r*(z). 

The  forward  and  backward  modes  a±(z,t)  can  be  obtained  by  (25).  From 
Vt±(z,t)  we  know  that  the  first  order  approximation  to  af(z,  <)  also  have  closed 
form  expressions. 

While  the  analysis  in  this  section  is  for  unit  step  excitation,  the  transient 
response  for  the  case  of  linear  terminal  conditions  to  general  input  waveforms  can 
be  easily  obtained  through  convolution. 


f 
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Figure  2.  Integration  paths  on  plane. 

SYSTEM  TERMINATED  WITH  NONLINEAR  LOADS 

Interconnecting  lines  in  integrated  circuits  are  usually  terminated  with  nonlinear 
loads,  e.g.,  the  input  impedances  of  transistors.  It  is  therefore  important  to  in¬ 
vestigate  the  transient  response  for  such  coupled  systems.  For  simplicity,  in  the 
following  analysis  we  assume  that  the  only  nonlinear  load  is  at  the  far  end  of  the 
active  line  (line  1  in  Fig.  1),  and  the  other  ends  of  the  coupled  line  system  are 
matched.  In  this  case,  the  boundary  conditions  become 


II 

o 

+  H 

a 

(40a) 

ar(^)  =  rn/(<)«J-(M) 

(406) 

e> 

p 

rt. 

II 

o 

(40c) 

o 

II 

d 

(40  d) 

We  also  assume  that  the  initial  states  for  all  independent  variables  are  zero. 

As  in  the  previous  section,  the  solutions  are  still  given  by  (24)  except  that  the 
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boundary  conditions  of  V  are 

-r«/(«i'(M)K w) 

ni  v*-n\*-y/i 

1 


vf{0  ,t)  = 


(<ll(0))2 

(<12(0)2 


‘r(M,+“VJ'(0’,) 

*2+(00  +  77T77TUr»<(ai'(O0)aJ'(f,<) 


(41a) 

(416) 

(41c) 

(41d) 


(<ll(0)2  2  '  ’  ;  (<ll(0)2 

We  will  not  elaborate  on  how  to  solve  V-*  from  the  nonlinear  equations  (41a)  - 
(41d) ,  but  wish  to  point  out  that  the  zeroth  order  solutions  to  (23)  will  consist 
not  only  of  the  terms  generated  by  the  excitation  but  also  the  reflections  due  to 
the  nonlinear  terminal  load.  The  expressions  of  the  zeroth  order  solutions  are 

^(xt*+(0))2-*1(0) 


v5.(0)(M>  -  2 


4- 


v/T?iMoyTMO)  +  v/iTMoJ^MO) 

n 

£  (S.(0)fl,(X)F  /((  -  r,(»)  -  2jT,) 

1=0 

ft(0)  £  [fi,(0)*i(<)P  %«  (\n-m) 


m= 0 


yrW(,  «>  _  2  {/(i^MQ))»-*»-(0) 

‘’n  V  '  v/1  ^  M°)  +  MO)  +  V1  ^*+(0)  -  MO) 

n 

53  {RiWXiWV  /(*  +  r*(*)  -  20  +  1)2<) 

7=0 

+  £  [«i(o)«i(oim 

m=0 

+  nW  -  2(m  + 1)15) 

(1  =  1,2) 


where 


rn*,«(z)  - 


(0  = 


1 

2 


2  ^/(lTMO)2  -fc-(0  rn/(z) 

y/1  =f  M*)  -  *-"(*)  +  v/iTM/r+  MO 

[(iuW)i^(0)«.')  -  (i.iwhn/V’') 
+(<uW)2*£(0,(<,<)  -  (<12(/))2Vj-(0)(<,1)] 


(  0,  if  t  <  2jTi  and  t  >  2 {j  +  1)1* 


(42) 

(43) 


(44) 


and  (<n)»  “d  (<12)*  ,  (*  =  1,2)  are  given  by  (29). 

Higher  order  terms  are  easier  to  calculate  because  a  linearized  model  can  be 
used  to  characterize  small  perturbations.  In  summary,  the  first  order  solutions 
have  the  following  form: 

In  region  ( nf )  :  2nT{  <  t £  <  (2n  +  l)Ti 

=  /.  dz>bi(z')VCn-\  (z'>*  +  T*(z)  ~  T»(z')) 

J*i.n 

4-  V~(1)  ( t+  1 

+  «,(«-l)6  V  «-n’  «,»/ 

(*U  *  ‘  +  n(<n))  (45a) 

*£/*(*.  0  =  J  dz'hi(z')Vi~n-\  (z'» *  “  rt(z)  +  Ti(z)) 

(2 nTi  <  t  <  2 tfn  -  2 nTt)  (456) 


In  region  ( n6 )  :  (2 n  +  1)7*  <  tin  <  2(n  +  1)TV 

+  lC/)  (v-^'n) 

(f-B<<<f-#  +  n(f)-r,(g)  (45c) 

VrW(z,  t)  =  j  dz'bi{z')vg0)  (z,  t  +  Tj(z)  -  Ti(z')) 

+  *(*)*£*  ^^  +  Ti(z)  -  T<)  +  ARnti 

((2 n  +  1)2*  <  *  <  2<~n  -  2nT,)  (45d) 

(<*1.2;  n  =  0,1,2,... ) 


where 


and 


A  +  Ar^(a^°»)«ti0)(<.'), 

<,  =  5  (*  +  +*(<  -  2nIi))  ■  ‘J.  “  |  K*  +  W)  +  n(z)| 

i  (»  -  r-‘(<  -  2(n  +  1)T,))  ,  l",  =  \  [(I  +  2(»  +  l)T.)  -  r,(z)] 
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The  remaining  steps  are  the  same  as  in  the  previous  section,  viz.  applying  (25) 
to  get 

•ft*.*)  =£((*ll(*))lVi+(*>t) 

“  (*  12(*))l  Vf  (z,  t)  ±  (tn(z))2V2- (z,  t)  T  {tl2(*)hV2  (z»  0]  (46a) 
«,_(M)  (*12(*))l  Vi4'(*i  *) 

+  (<ll(*))l  V- (z,  t)  T  (tU{z))2V+(z,  t)  ±  (tn(2))2V{(z, «))  (466) 

(i  =  l,2) 

NUMERICAL  RESULTS 

To  verify  our  formulation,  we  first  examine  &  special  case  where  fee  =  •  The 

parameters  of  Fig.  3  are  taken  from  Fig.  14  of  [19],  i.e.  —  *(*)  = 

0.45  +  0.2tanh(2.5z  —  1.25).  Both  Zt  and  Z 1  are  matched  to  the  line,  and 
the  units  of  distance  and  time  are  normalized  to  the  lengths  of  the  lines  and  the 
corresponding  signal  travel  times,  respectively.  The  excitation  es(t)  is  assumed 
to  be  a  pulse  of  amplitude  4  with  rise  time  of  0.25T,  fall  time  of  O.lT,  and 
duration  of  0.15T  *  The  first  order  approximation  shown  here  matches  very  well 
with  the  result  in  [19]. 

The  power  of  the  present  approach  is  in  its  capability  of  computing  codirec- 
tional  as  well  as  contra-directional  coupling.  For  the  rest  of  the  calculation, 
we  assume  that  the  transmission  line  is  in  a  medium  with  relative  permittiv¬ 
ity  «r  =  9  and  line  length  10  mm  and  the  input  signal  is  a  pulse  of  amplitude 
5.  Figures  4  through  7  show  the  far-end  crosstalks  when  both  lines  have  linear, 
matched  loads.  For  Figs.  4  and  5,  kc{z)  =  0.15  +  0.05sin(0.1irz  +  w/4)  and 
k i{z )  =  0.10  +  0.05sin(0.l7rz  4-  w/4)  where  the  unit  of  length  is  mm.  Hence  the 
co- directional  coupling  coefficient  i fe+  =  (k(j  —  ki)/ 2  =  0.025  is  a  constant.  The 
only  difference  is  that  the  rise  time  and  the  fall  time  are  10  ps  for  Fig.  4,  but  20 
ps  for  Fig.  5.  The  durations  of  the  pulses  are  also  different:  20  ps  for  Fig.  4  and 
30  ps  for  Fig.  5.  The  results  show  good  agreement  with  the  approximate  formula 
in  [12].  Indeed,  the  far-end  crosstalk  is  inversely  proportional  to  the  rise  time 
since  k( 7  and  ki  are  so  close  that  the  difference  between  end-to-end  traveling 
times  of  modes  Vf  and  V2  (or  Vj”  and  V^-  )  *s  ^e8a  tban  the  rise  time.  Note, 
however,  that  the  formula  in  [12]  is  rendered  invalid  when  fc+(z)  is  no  longer 
uniform.  Figures  6  and  7  both  have  ki  —  0.25  +  0.10tanh(0.25z  —  1.25).  k q 
is  0.25  -I-  0.15tanh(0.25z  —  1.25)  for  Fig.  6  and  0.25  +  0.23tanh(0.25z  -  1.25) 
for  Fig.  7.  Equivalently,  t+(z)  =  0.025  tanh(0.25z  —  1.25)  for  Fig.  6  and  k+  = 
0.05  tanh(0.25z  -  1.25)  for  Fig.  7.  Qualitatively,  the  results  are  consistent  with 
the  formula  in  [12j.  However,  the  peak  values  are  significantly  smaller  than  those 
predicted  by  the  previous  approximate  formula  if  the  maximum  k+(z)  is  sub¬ 
stituted.  This  is  understandable,  since  the  positive  and  negative  parts  of  fc+(z) 
tend  to  average  out. 

*  The  rise  time  and  the  fall  time  are  mistakenly  written  as  0.125T  and  0.0625T 
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In  the  following  figures  we  consider  systems  having  coupling  coefficients  of  the 
form: 

A  k 

k-(z)  =  k0  +  —[1  -  cos (az  +  /3)j  (47a) 

it 

and 

A  i 

M*)  =  ~ 2~[1  -  cos(otz  +  /?)]  (47 b) 

and  a  nonlinear  load  at  the  far-end  of  the  active  line  (port  3  in  Fig.  1)  characterized 
by 

rnlia\ )  =  rm[l  +  tanh(7(o!  -  a0))]  +  r0  (48) 

where  a\  =  -\-a^(£,t).  The  remaining  ports  are  assumed  to  be  matched. 

The  input  signal,  as  shown  in  Fig.  8,  has  both  rise  time  and  fall  time  of  10  ps, 
duration  20  ps  and  amplitude  5.  The  parameters  of  Fig.  8  are  given  as  follows: 

k0  =  0.2,  A k  =  —0.05,  a  =  2ir  rad /mm,  /3  =  0 

and 

rm  —  -0.3,  7  =  1.0,  ao  =  2.5  and  r0  =  -rm(l  —  tanh(7a0)) 

In  Fig.  8,  the  first  order  transmitted  signal,  reflected  signal,  near-end  and 
far-end  crosstalks  are  displayed  together  to  allow  comparison  of  their  relative 
magnitudes.  Figs.  9-12  show  the  differences  between  the  zeroth  order  and  the 
first  order  solutions  for  the  reflected,  transmitted  signal,  near-end  and  far-end 
crosstalks.  The  small  discrepancies  justify  the  use  of  the  iterative  method. 

It  is  of  interest  to  see  what  happens  when  Je+(z)  is  positive  instead  of  negative. 
In  Figs.  13-17,  we  choose  A k  =  +0.05,  and  the  coupling  coefficients  become 

fc_(z)  =  0.2  +  0.025[1  -  cos(27rz)] 
k+(z)  =  0.025(1  —  cos(27rz)] 

The  responses  at  all  four  ports  are  shown  in  Fig.  13  whereas  Figs.  14-17  provide 
comparison  with  the  previous  example  in  which  A k  =  —0.05.  We  notice  that 
the  far-end  crosstalk  for  this  case  is  nearly  the  negative  of  that  for  the  previous 
case.  This  can  be  explained  in  terms  of  the  differences  in  mode  velocities.  From 
the  theory  of  coupled  transmission  lines  one  can  show  that  the  far-end  crosstalk 
is  equal  to  the  difference  of  odd  and  even  modes,  which  correspond  to  our  V* 
and  V ^  .  The  velocity  of  the  even  mode  may  be  larger  or  smaller  than  that  of 
the  odd  mode,  depending  on  whether  is  positive  or  negative.  Therefore,  the 
initial  arrival  of  the  faster  mode  at  the  far  end  results  in  positive  crosstalk  in  one 
case  and  negative  in  the  other.  The  deviation  after  the  first  arrival  (  ~  lOOps )  is 
due  to  the  nonlinear  load. 

CONCLUSIONS 

General  approximate  solutions  to  the  transient  response  on  two  identical,  nonuni- 
formly  coupled  transmission  lines  terminated  with  linear  or  nonlinear  loads  have 
been  obtained  through  an  iterative  scheme  that  was  previously  applied  to  a  smaller 
class  of  problems  in  [19].  The  iterative  method  is  very  useful  when  the  coupling 


Figure  5. 


Response  of  line  2  if  kc(z)  =  0.15  +  0.05  sin(0.1irz  +  ir/4)  and 
kl(z)  =  0.10  +  0.05sin(0.l7T2  +  ir/4) .  tr  —  20  ps 


Figure  6.  Response  of  line  2  if  kc(z)  =  0.25  +  0.15  tanh(0.25r  -  1.25)  and 
kL{z)  =  0.25  +  0.10  tanh(0.25r  -  1.25).  tr  =  10  ps 
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Figure  7. 


Figure  8. 


Time  i  p« ) 


Response  of  line  2  if  k(j(z)  =  0.25 +  0.20  tanh(0.25z  -  1.25)  and 
k l(z )  =  0.25  +  0.10  tanh(0.25z  -  1.25) .  tr  =  10  ps 
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Response  at  different  ports  of  a  system  with  k-(z)  —  0.20  — 
0.025(1  —  cos27rz)  and  &+(z)  =  —0.025(1  -  cos2wz). 
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Figure  9.  Comparison  of  the  zeroth  and  the  first  order  solutions  of  the 
reflected  signal,  i fc_(z)  =  0.2  -  0.025(1  -  cos2irz),  *+(z)  = 
-0.025(1  -  cos2ttz) 


Figure  10.  Comparison  of  the  zeroth  and  the  first  order  solutions  of  the 
transmitted  signal.  k-(z)  =  0.2  -  0.025(1  —  cos2irx) ,  fc+(z)  = 
-0.025(1  -  cos2ttz) 
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Figure  11.  Comparison  of  the  zeroth  and  the  first  order  solutions  of  the 
near-end  crosstalk.  k-(z)  =  0.2  -  0.025(1  -  cos  2 ttz)  ,  fc+(z)  = 
-0.025(1  —  cos2irz) 


TIMEC  PS) 

Figure  12.  Comparison  of  the  zeroth  and  the  first  order  solutions  of  the 
far  end  crosstalk.  k-(z)  =  0.2  -  0.025(1  -  cos2ttz)  ,  k+(z)  = 
-0.025(1  —  cos2rrz) 


/ 
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Figure  13.  Response  at  different  ports  of  a  system  with  Jfe_(r)  =  0.20  + 
0.025(1  -  cos  27rr)  and  &+(z)  =  +0.025(1  -  cos  2m). 


Figure  14.  Comparison  of  reflected  signals  on  lines  terminated  nonlinear 
loads  with  fc+(z)  of  opposite  signs. 
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Figure  15.  Comparison  of  transmitted  signals  on  lines  terminated  nonlin¬ 
ear  loads  with  ifc+(z)  of  opposite  signs. 
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Figure  16.  Comparison  of  near  end  crosstalk  on  lines  terminated  nonlinear 
loads  with  lfe+(z)  of  opposite  signs. 
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Figure  17.  Comparison  of  far  end  crosstalk  on  lines  terminated  nonlinear 
loads  with  fc+(z)  of  opposite  signs. 

coefficients  are  slowly  varying  with  position  since  the  zeroth  order  or  first  order 
approximation  would  be  sufficiently  accurate  yet  much  easier  to  calculate.  Fur¬ 
thermore,  with  the  help  of  newly  devised  special  transformations,  we  have  shown 
that  both  the  codirectional  coupling  and  contradirectional  coupling  of  the  problem 
with  unit-step  excitation  and  linear  loads  have  closed  form  expressions  up  to  the 
first  order  approximation.  Arbitrary  excitation  can  then  be  taken  into  account 
by  convolution.  This  method  is  hence  most  efficient.  As  for  nonlinear  termina¬ 
tions,  numerical  integrations  are  performed  along  the  characteristics.  Examples 
have  been  given  for  both  cases  to  illustrate  the  use  of  this  method.  Extension  to 
problems  in  which  the  phase  velocities  of  coupled  lines  are  not  equal,  or  where 
more  than  two  coupled  lines  are  involved  is  also  under  consideration. 
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